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This  study  examines  the  use  of  genetic  search  techniques  in  the 
multicriterion  optimal  design  of  structural  systems.  The  global  nature  of  these 
search  techniques  makes  them  better  suited  for  design  spaces  that  are  known  to 
be  disjoint  or  nonconvex.  The  study  was  motivated  by  an  application  involving  the 
optimal  design  of  viscoelastically  damped  composite  structures.  In  addition  to  the 
multicriterion  nature  of  this  design  problem,  its  design  space  is  a  mix  of 
continuous,  integer,  and  discrete  design  variables.  The  more  traditional 
mathematical  programming  based  methods  of  optimization  are  known  to  generate 
suboptimal  designs  for  this  class  of  problems. 

Genetic  algorithms  belong  to  a  general  category  of  stochastic  search 
methods.  An  optimal  design  is  located  by  combining  favorable  characteristics  of 
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several  designs  in  such  a  way  as  to  maximize  a  defined  fitness  function.  The 
transformation  operators  that  facilitate  this  combination  have  their  philosophical 
basis  in  biological  principles  of  evolution.  As  shown  in  this  work,  the  method  can 
easily  accommodate  the  presence  of  integer  and  discrete  variables  in  the  design 
space.  An  extension  of  genetic  search  to  the  multicriterion  optimal  design 
problem  is  presented  in  this  work.  This  design  problem  requires  that  an  optimal 
design  be  generated  for  which  each  of  the  criteria  are  simultaneously  optimized. 
Two  distinct  strategies  for  generating  such  designs  are  presented  in  this  work.  A 
distinct  feature  of  these  strategies  is  that  they  allow  for  a  simultaneous  generation 
of  optimal  designs  in  problems  where  the  individual  criterion  are  weighted 
differently. 

The  design  of  viscoelastically  damped  composite  beam  structures 
considered  in  this  work  was  based  on  a  three-dimensional  finite  element  solution. 
This  analysis  is  computationally  demanding,  and  in  order  for  it  to  be  used  for 
optimization,  approximations  to  the  exact  analysis  were  required.  A  multilayer 
perceptron  model  of  an  artificial  neural  network  was  used  to  map  the  input- 
output  relationships  between  the  design  variables  and  the  output  response,  and 
this  approximation  was  used  in  the  optimal  design  in  lieu  of  the  exact  analysis. 
Optima]  designs  of  composite  beams  for  minimum  weight  and  maximum  damping 
were  obtained  using  this  approach. 
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CHAPTER  1 
INTRODUCTION 

1.1  Literature  Survey 

The  earliest  theoretical  developments  in  the  field  of  structural  optimization 
may  be  traced  back  to  a  publication  by  Maxwell  [1]  in  1869.  In  1904,  Michell  [2] 
proposed  a  theory  which  provided  valuable  insights  into  determining  the  optimum 
configuration  of  truss  structures  under  a  single  load;  however,  the  approach  had 
limited  practical  applications.  In  the  1940's  there  was  a  new  development  in  the 
theory  for  optimal  design  of  structural  components--the  simultaneous  failure  mode 
theory.  As  the  name  suggests,  the  structural  components  were  sized  in  a  manner 
such  that  failure  would  occur  simultaneously  in  several  failure  modes.  In  the  area 
of  plastic  collapse  design,  linear  programming  methodology  was  used  in  the  early 
1950's  to  obtain  the  optimum  configuration  of  truss  and  frame  structures. 

Since  1960,  optimality  criteria  methods  and  nonlinear  programming 
algorithms  (NLP)  have  become  two  of  the  most  popular  structural  optimization 
methods.  The  optimality  criterion  methods  belong  to  a  class  of  indirect  methods 
for  optimal  design.  A  structure  with  given  loading  conditions  and  prescribed 
constraints  is  deemed  optimal  when  its  derived  optimality  criterion  is  satisfied  via 
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numerical  algorithms.  This  optimality  criterion  technique  [3-5]  is  not  general,  as 
the  optimality  criterion  definition  is  problem  dependent,  and  quite  often  difficult 
to  obtain.  The  use  of  nonlinear  programming  techniques  in  optimum  structural 
design  affords  a  more  direct  and  flexible  solution  procedure.  Professor  L.  A. 
Schmit  pioneered  this  approach,  now  being  widely  recognized  as  a  viable  design 
tool  in  the  industry.  A  number  of  key  review  publications  and  illustrations  of 
applications  are  available  in  References  6-17. 

Although  NLP  methods  have  been  developed  to  the  point  where  they  can 
be  used  almost  routinely,  they  have  several  drawbacks.  Such  drawbacks  may  result 
from  having  an  optimization  problem  with  a  large  number  of  design  variables  and 
constraints.  Furthermore,  the  design  space  may  be  disjoint  or  nonconvex. 
Nonlinear  programming  methods  are  not  suitable  for  an  optimization  problem 
involving  one  or  more  design  variables  which  are  of  discrete/integer  type  because 
of  the  difficulties  in  applying  the  one-dimensional  search  technique,  and  in  the 
determination  of  gradients  of  objective  function  and  constraints.  A  combinatorial 
programming  strategy  has  been  proposed  to  overcome  the  discreteness  of  the 
design  space.  This  technique  demands  significant  investment  of  computational 
resources  which  worsens  as  the  dimensionality  of  the  non-continuous  design 
variables  increases.  Since  the  traditional  NLP  based  methods  essentially  represent 
a  local  search  strategy,  they  can  only  determine  the  relative  optimum  located 
nearest  to  the  starting  design.  Further,  as  the  more  efficient  methods  belonging  to 
the  category  of  NLP  require  at  least  the  first  order  gradient  information  of  the 
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objective  function  and  constraints,  they  will  very  likely  fail  as  the  search  passes  the 
disconnected  region  of  a  disjoint  design  space. 

Decisions  in  engineering  design  typically  involve  multiple  and  frequently 
conflicting  requirements.  Perhaps  the  first  effort  in  the  field  of  multicriterion 
optimization  may  be  attributed  to  an  Italian  economist,  Pareto,  who  in  1896 
introduced  the  concept  within  the  framework  of  welfare  economics  [18].  The 
implications  of  his  work  to  optimization  theory,  operations  research,  and  control 
theory  were  only  recognized  in  the  late  1960's.  Multicriterion  optimization  in 
engineering  design  has  been  examined,  and  numerous  examples  of  the  design  of 
mechanical  structures  have  been  presented  [19-27].  Despite  the  multicriterion 
nature  of  most  design  problems,  the  principal  focus  of  most  research  efforts  has 
been  on  pursuing  an  efficient  algorithm  for  the  scalar  objective  function 
optimization  problem.  One  common  approach  to  solving  multicriterion 
optimization  problems  is  to  select  one  of  the  criteria  as  the  objective  function,  and 
to  trade  off  among  the  remaining  criteria  using  appropriate  design  constraints. 
The  apparent  simplicity  afforded  by  this  method  is  attractive;  however,  it  should 
be  noted  that  the  treatment  of  criteria  as  constraints  does  not  yield  the  same 
optimum  design  as  would  be  obtained  when  solving  the  optimum  design  problem 
as  one  possessing  multiple  objectives  [28]. 

The  goal  of  this  research  was  to  develop  efficient  methods  for 
multicriterion  optimal  design.  The  approach  was  then  applied  to  optimally  design 
viscoelastically  damped  structures  with  maximum  damping  and  minimum  weight. 
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This  exactly  characterizes  a  multicriterion  structural  optimization  problem.  In  this 
problem,  however,  the  multicriterion  optimization  problem  was  complicated  by 
the  presence  of  discontinuity  in  the  design  space.  The  design  variables  were  a  mix 
of  integer,  discrete,  and  continuous  types.  Integer  and  discrete  design  variables 
were  necessary,  because  the  number  of  composite  laminae  forming  a  laminate  is 
characterized  by  an  integer  variable,  and  the  thicknesses  of  the  constraining  and 
viscoelastic  layers  were  restricted  to  a  manufacturer-provided  list  of  discrete 
values.  This  multicriterion  optimal  design  problem  comprises  a  complicated 
design  space  which  is  discontinuous,  and  nonconvex. 

A  finite  element  based  analysis  of  the  structural  damping  ability  of  a 
constrained-layer  treated  composite  beam  was  investigated  by  Sun  and  his  co- 
workers [29-33].  The  loss  factor  of  a  composite  laminate  of  arbitrary  stacking- 
sequence,  and  with  a  discrete  type  of  add-on  viscoelastic  damping  treatment  can 
be  effectively  and  accurately  predicted.  This  work  was  used  in  the  present 
research  as  an  analysis  tool  for  determining  the  characteristics  and  the 
dynamic/static  response  of  viscoelastically  damped  composite  structures. 

In  order  to  solve  the  multicriterion  optimization  problem  with  a 
discontinuous  and/or  nonconvex  design  space  effectively,  two  methods  were 
explored.  The  first  approach  was  an  NLP  based  procedure  with  a  branch-and- 
bound  type  algorithm  used  to  account  for  the  integer  and  discrete  design 
variables.  The  second  approach  was  the  use  of  genetic  algorithms,  and  this 
comprises  a  principal  focus  of  the  present  study. 
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Genetic  algorithms  were  first  proposed  by  John  Holland  in  1975  [34],  and 
have  since  been  adapted  for  a  large  number  of  applications  in  game  theory, 
induction  systems,  and  other  aspects  of  human  cognition  such  as  pattern 
recognition  and  natural  language  processing  [35-38].  These  algorithms  are 
generally  regarded  to  be  in  the  same  category  of  stochastic  search  methods  as 
simulated  annealing  [39].  Both  approaches  have  their  basis  in  natural  processes 
(simulated  annealing  is  derived  from  the  principles  of  statistical  mechanics). 

Genetic  search  methods  are  based  on  Darwin's  theory  of  survival  of  the 
fittest  [40].  A  set  of  design  alternatives,  (analogous  to  a  population  in  a  given 
generation)  are  allowed  to  reproduce  and  cross  among  themselves,  with  bias 
allocated  to  the  most  fit  members  of  the  population.  Combinations  of  the  most 
desirable  characteristics  of  the  mating  pairs  of  the  population  result  in  progenies 
that  are  more  fit  than  either  of  the  parents.  If  a  measure  which  indicates  the 
fitness  of  a  generation  is  also  the  desired  goal  of  a  design  process,  successive 
generations  produce  better  values  of  the  objective  function.  An  obvious 
advantage  in  this  approach  is  that  the  search  is  not  based  on  gradient  information, 
and  has  no  requirements  on  the  continuity  or  convexity  of  the  design  space. 
Another  advantage  is  that  the  solution  space  in  genetic  algorithms  is  a  discrete 
mapping  of  the  design  space;  thereby  a  design  space  with  a  mk  of  integer  and 
discrete  design  variables  can  be  treated  efficiently. 

While  the  potential  of  genetic  algorithms  as  a  function  optimizer  has  been 
a  subject  of  more  recent  interest,  the  first  use  of  genetic  search  methods  in  the 


context  of  structural  optimization  was  done  by  Hajela  [41,42].  In  that  work  [42], 
globally  optimal  structural  designs  were  generated  using  genetic  search  methods, 
despite  the  existence  of  a  nonconvex  and  disjoint  design  space. 

The  present  research  proposes  and  evaluates  a  number  of  distinct  coding 
strategies  for  accommodating  integer,  discrete,  and  continuous  design  variables  in 
genetic  algorithms.  It  also  extends  the  genetic  search  method  to  simultaneously 
obtain  multicriterion  optimal  designs  for  several  cases  where  the  individual 
criterion  are  weighted  differently  relative  to  the  remaining  criteria. 

1.2  Scope  of  Present  Research 

The  scope  of  the  present  work  and  its  organization  can  be  summarized  as 
follows.  Damping  in  composite  structural  elements  is  the  subject  of  chapter  two. 
Both  the  inherent  damping  characteristics  of  composite  structures  and  add-on 
types  of  passive  damping  are  discussed.  The  influences  of  constituent  material 
properties  and  the  geometrical  configuration  of  the  composite,  on  its  damping 
characteristics,  are  presented.  Three  optimal  design  problems  for  damping 
enhancement  are  formulated,  and  a  detailed  implementation  of  these  design 
problems  is  presented. 

Chapter  three  provides  an  introduction  to  the  traditional  nonlinear 
programming  techniques,  with  a  focus  on  the  feasible-usable  directions  method. 
This  method  was  used  in  conjunction  with  a  modified  branch-and-bound  strategy 
for  solving  discrete/integer  design  variable  problems.  The  min-max  global 
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criterion  strategy  for  multicriteria  optimization  is  introduced.  A  two-dimensional 
idealization  of  a  cantilever  beam  with  add-on  viscoelastic  damping  layers,  is 
optimally  configured  as  an  example  of  an  NLP  based  method  for  multicriterion 
optimal  design  of  structural  sandwich  beam. 

In  chapter  four,  the  theoretical  background  and  the  basic  elements  of 
genetic  search  methods  are  reviewed.  The  binary  implementation  of  integer, 
discrete,  and  continuous  design  variables  is  discussed.  Two  distinct  approaches  for 
generating  multicriterion  optimal  designs  are  presented  in  this  chapter.  These 
methods  are  proposed  as  an  alternative  to  the  more  computationally  intensive 
NLP  based  schemes  for  multicriterion  design. 

Chapter  five  describes  the  numerical  implementation  of  genetic  search 
methods  described  in  chapter  four.  To  demonstrate  the  use  of  genetic  algorithms 
in  structural/mechanical  optimization  problems  with  a  mix  of  integer,  discrete, 
and  continuous  design  variables,  a  lap  rivet  efficiency  maximization  problem,  a 
Belleville  spring  sizing  and  stacking  problem,  and  an  optimal  structural  sizing 
problem  were  investigated.  The  optimal  stacking-sequence  of  square  and 
rectangular,  simply-supported  laminated  composite  plates,  for  maximum  transverse 
stiffness,  was  also  considered.  Genetic  search  methods  for  simuhaneously 
generating  multicriterion  optimal  designs  were  investigated  on  a  number  of 
structural  optimization  problems  with  multiple  design  criteria.  Included  in  this 
study  was  a  viscoelastically  damped,  laminated  composite  beam  problem  which 
used  maximum  damping  and  minimum  weight  as  optimization  criteria. 
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In  chapter  six,  a  neural-network-based  approximation  of  structural  analysis 
is  proposed  for  use  in  conjunction  with  genetic  algorithms.  This  significantly 
alleviates  the  required  computational  effort  in  the  optimal  design.  A  number  of 
multilayer  perceptron  models  of  neural  networks  and  the  back-propagation 
algorithm  used  in  the  training  of  such  networks,  are  also  discussed.  This 
approximation  strategy  is  examined  in  the  context  of  a  multicriterion  optimization 
problem  where  a  3-dimensional  cantilever  composite  beam  with  add-on 
viscoelastic  damping  tape  is  optimally  configured  for  maximum  damping.  The 
influence  of  such  damping  mechanisms  on  the  viscoelastically  damped  cantilever 
beam  is  studied  in  this  implementation. 

Conclusions  drawn  from  this  study  and  recommendations  for  future 
research  are  presented  in  the  final  chapter. 


CHAPTER  2 
DAMPING  TREATMENT  AND  ANALYSIS 

2.1  Background 

The  need  for  efficient  light-weight  structures  in  the  aerospace  industry  has 
dramatically  increased  the  use  of  composite  materials  in  the  last  decade.  The 
application  of  composite  materials  in  structures  that  are  subjected  to  dynamic 
loads  has  introduced  stringent  requirements  on  the  dynamic  characteristics  of 
these  materials.  Vibration-induced  damage  in  structures  has  long  been  one  of  the 
leading  causes  for  structural  failure.  Scientists  and  engineers  have  dedicated 
tremendous  efforts  to  controlling  and  reducing  vibration. 

Active  and  passive  control  are  among  the  most  popular  approaches 
currently  used  to  improve  the  damping  ability  of  structures.  Active  control  utilizes 
external  devices  referred  to  as  sensors  and  actuators  to  shape  the  dynamic 
characteristics  of  the  system.  On  the  other  hand,  passive  control  mechanisms 
emphasize  the  integrated  optimal  design  of  inherent  structural  characteristics  and 
external  damping  devices  to  maximize  damping  under  expected  dynamic  loads. 
More  specifically,  for  a  laminated  structural  composite,  an  optimal  passive 
damping  can  be  achieved  by  a  combination  of  configuring  the  layup  and  material 
selection  of  the  base  structure  (the  laminated  composite  itself),  and  by  the 
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addition  of  external  damping  materials  such  as  viscoelastic  tapes  to  the  surface  of 
the  laminated  composite.  The  optimization  problem  associated  with  tailoring  the 
passive  damping  ability  of  viscoelastically  damped  composite  structural  elements 
was  the  main  focus  of  this  study. 

Damping  is  the  ability  of  a  structure  to  suppress  the  vibration  through  the 
dissipation  of  vibratory  energy.  In  earth-bound  structures,  there  are  a  number  of 
damping  mechanisms  in  which  vibratory  energy  is  converted  to  other  forms  of 
energy  such  as  sound  or  heat  [43].  These  mechanisms  are  classified  under  the 
class  of  non-material  damping.  Energy  can  also  be  dissipated  as  the  material 
deforms,  and  in  such  cases  the  damping  mechanism  is  active  within  the  material 
volume  [43-45].  Material  damping  is  also  referred  to  as  either  internal  damping  or 
structural  damping.  This  damping  is  caused  by  the  energy  dissipated  due  to 
mechanisms  such  as  viscoelastic  behavior  in  polymeric  material  or  the  interfacial 
slip  in  metallic  materials. 

The  principal  mechanism  for  damping  in  glass  and  graphite  reinforced 
polymer  matrk  composites  is  the  viscoelastic  energy  dissipation  that  occurs  in  the 
polymer  matrix.  Under  a  normal  stress  level,  the  composite  exhibits  linear 
viscoelastic  behavior,  and  its  dynamic  stiffness  and  damping  are  independent  of 
vibration  amplitude. 

The  first  attempt  to  tailor  the  damping  ability  of  basic  composite  structures 
can  be  traced  to  an  effort  due  to  Plunkett  and  Lee  [46].  Several  models  have  been 
proposed  for  predicting  the  damping  ability  in  short  fiber  reinforced  composites. 
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Sun  and  his  co-workers  [47]  have  developed  theoretical  relationships  for  the 
material  damping  of  aligned  short  fiber-reinforced  polymer  matrix  composites 
under  off-axis  loading.  The  complex  modulus,  E^',  of  the  short-fiber  composite 
can  be  expressed  as, 

1  ^       1       ^  cos'^e  ^  sin^e 

^  '^x 

where  E^'  and      are  the  storage  and  loss  moduli,  respectively;  E*l,  E'-j-  and  G'l- 
represents  the  complex  form  of  the  elastic  modulus  in  the  longitudinal  and 
transverse  directions,  and  the  shear  modulus,  respectively;  v^-j-  represents  the  in- 
plane  Poisson's  ratio.  The  damping  or  loss  factor  can  be  obtained  as  follows 


(2.2) 


sin^0cos-0 


(2.1) 


The  model  based  on  Cox's  shear  lag  theory,  as  shown  in  Figure  2.1,  allows 
us  to  express  the  loss  factor  (a  measure  of  the  structural  damping)  as  a  function 
of  fiber  and  matrix  properties  and  the  geometric  parameters  shown  in  Figure  2.1. 

E" 

=  vA^n  El  Gl  p,  s,  d,  D,  6)  (2.3) 


where  E'f,  E\^,  and  G'„,  represent  the  complex  form  of  the  modulus  of  the  short 
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Figure  2.1  Representative  volume  element  for  off-axis  loading. 
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fiber,  the  modulus  of  matrix,  and  the  shear  modulus  of  matrix,  respectively.  A 
detailed  expression  for  this  is  presented  by  Gibson  et  al.  [48].  Experiments  [49] 
indicate  that  the  damping  ability  in  composites  reinforced  by  discontinuous  fibers 
is  generally  greater  than  that  in  continuous  fiber  composites.  Sun  and  Wu  [50] 
attributed  this  to  the  stress  concentration  effects  present  at  the  fiber  ends,  which 
facilitate  the  transfer  and  dissipation  of  energy  in  the  viscoelastic  polymer  matrix. 
The  analysis  [48,50]  clearly  indicated  the  dependence  of  internal  material  damping 
in  composites  on  the  fiber  and  matrix  geometry.  Parameters  such  as  loading  angle, 
fiber  aspect  ratio,  volume  fraction,  fiber  length,  fiber  diameter,  etc.  were  all  found 
to  influence  the  dynamic  performance  of  the  composite.  The  analysis  also  showed 
that  the  increase  in  damping  ability  was  usually  associated  with  a  corresponding 
stiffness  loss.  A  parametric  study  can  be  conducted  by  varying  each  of  these 
variables,  one  at  a  time,  to  examine  their  influence  on  the  damping.  However, 
parametric  studies  of  this  form  have  been  shown  to  yield  suboptimal  results. 
Hajela  and  Shih  [51]  proposed  a  design  synthesis  approach  based  on  formal 
multicriterion  optimization.  Two  distinct  analysis  models.  Cox's  shear-lag  theory 
and  an  advanced  shear-lag  model,  were  used  in  their  work.  With  this  approach, 
the  extensional  loss  factor  of  a  representative  volume  element  (RVE)  was 
maximized  subject  to  constraints  on  the  element  mass  and  stiffness  characteristics. 

Damping  in  laminated  fiber  composites  can  be  determined  through  the 
laminated  plate  theory  approach.  The  in-plane  material  damping,  r]-^,  of  the 
laminated  composite  is  defined  as  follows. 
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E" 

ri,  =  ^  ij  =  1,2,6  (2.4) 

F 

where  E",j,  and  E'jj  represent  the  real  and  imaginary  part  of  the  in-plane  complex 
moduli  E'ij,  which  can  be  expressed  as 

^(z  =  —r^,  =  ^.j  ^  ^^iJ  (2.5) 

(4) 

where  (A  -j)"^  is  the  element  in  the  i'"  row  and  j"'  column  of  the  inverse  of  matrix 
[A'ij],  which  is  the  complex  form  of  the  reduced  in-plane  stiffness  matrix,  [Ajj]. 
Damping  in  the  laminated  fiber  composites  can  also  be  determined  through  an 
energy  approach  in  conjunction  with  a  3-dimensional  finite  element  method.  This 
model  represents  a  more  realistic  approach  by  including  the  energy  dissipated  at 
the  interfaces.  Detailed  derivation  of  damping  with  both  approaches  is  available 
in  [52]. 

Although  the  internal  damping  capabiHty  in  laminated  composites  can  be 
adjusted,  the  extent  to  which  it  can  be  increased  is  limited  by  other  design 
requirements.  To  further  increase  the  damping  ability,  viscoelastic  surface 
damping  treatments  can  be  used.  Add-on  viscoelastic  treatments  can  be  used  to 
increase  the  damping  ability  of  the  laminated  composites  without  decreasing  the 
stiffness  of  the  laminated  composites.  Extensional  damping  treatment,  a  commonly 
used  surface  damping  treatment,  is  also  referred  to  as  unconstrained-  or  free-layer 
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damping  treatment.  The  viscoelastic  material  is  coated  on  one  or  both  sides  of  a 
structure;  whenever  the  structure  is  subjected  to  cyclic  bending,  the  damping 
material  will  be  subjected  to  tension-compression  deformation.  The  synchronized 
tension-compression  deformation  of  the  viscoelastic  layer  dissipates  the  energy 
away  from  the  vibrating  base  structure. 

The  shear  type  of  damping  treatment  provides  a  more  efficient  damping 
enhancement  ability  than  the  unconstrained  layer  treatment  on  a  constant  weight 
basis  [53].  This  is  largely  due  to  the  fact  that  the  energy  dissipation  in  the 
viscoelastic  material  is  entirely  due  to  the  shear  deformation  and  is  independent 
of  dilatational  deformation  [54].  The  shear  damping  treatment  is  similar  to  the 
unconstrained-layer  type,  except  that  the  viscoelastic  material  is  constrained  by  a 
stiff  constraining  layer.  Therefore,  whenever  the  structure  is  subjected  to  cyclic 
bending,  the  constraining  layer  will  force  the  viscoelastic  material  to  deform  in 
shear.  This  shear  deformation  is  the  mechanism  by  which  energy  is  dissipated. 
Figure  2.2  shows  a  sketch  of  both  free  and  constrained  layer  damping  treatments. 
It  has  been  noted  that  the  overall  mass,  stiffness,  and  damping  characteristics  of  a 
structure  will  determine  its  damping  abihty.  All  of  the  above  three  parameters, 
and  not  just  the  damping  material,  need  to  be  optimally  configured  in  order  to 
obtained  the  highest  damping  ability  under  some  design  constraints.  Typically,  in 
laminated  composite  structures,  several  factors,  such  as  the  stacking  sequence  in 
the  laminate,  the  location,  amount,  and  type  of  treatment,  strongly  influence  the 
damping  response. 

i 
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Figure  2.2  Illustration  of  mechanisms  of  damping  treatments. 


17 


2.2  Measures  of  Damping 

There  are  various  damping  measures  which  are  commonly  used  in  different 
areas.  Three  common  definitions  of  damping  are  introduced  in  this  section. 

For  a  single  degree  of  freedom  system  shown  in  Figure  2.3,  with  viscous 
damping  and  an  external  excitation  F(t),  the  differential  equation  of  motion  is, 


where  m  represents  the  mass  of  the  system,  c  represents  the  viscous  damping,  and 
k  represents  the  spring  constant.  For  viscously-damped  free  vibration,  the  general 
solution  to  the  homogeneous  equation  is 


mx  +  cx  +  kx  =  F{t) 


(2.6) 


(2.7) 


where 


11 


a  = 


c 


k 


(2.8) 


m 


Cj  and  c,  are  constants  to  be  determined  by  initial  conditions.  An  important 
characteristic  of  a  structural  system,  the  damping  ratio  or  damping  factor,  is 
defined  as 
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Figure  2.3  Sketch  of  viscous  damping  model. 
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where 

=  2  sf^  (2.10) 

A  value  c^,  for  c  determines  the  system  to  be  overdamped  (c  >  cj,  underdamped 
(c  <  c^),  or  critically  damped  (c  =  cj.  Depending  whether  it  is  overdamped, 
imderdamped,  or  critically  damped,  the  system  exhibits  an  oscillatory  motion, 
nonoscillatory  motion,  or  critically  damped  motion,  respectively. 

A  convenient  way  to  determine  the  amount  of  damping  present  in  a  system 
is  to  measure  the  rate  of  decay  of  free  oscillation.  The  logarithmic  decrement,  6, 
is  defined  as  the  natural  logarithm  of  the  ratio  of  any  two  successive  amplitudes  of 
a  free  vibration,  as  shown  in  Figure  2.4.  A  relation  exists  between  the  damping 
ratio,  I,  and  the  logarithmic  decrement,  <S,  under  the  assumption  of  low  damping; 
the  relation  is  $  =  8 /it. 

A  commonly  used  definition  of  structural  damping,  also  referred  to  as 
hysteretic  damping  was  originally  formulated  as  a  basis  for  describing  the  internal 
damping  properties  of  solids  [55].  If  the  hysteresis  loop  for  a  material  has  a 
certain  shape,  the  damping  can  be  considered  by  using  a  complex  stiffness,  K'  ( = 
K'  +  iK"),  where  the  real  and  imaginary  part  are  referred  to  as  the  storage 
modulus  and  loss  modulus,  respectively.  The  loss  factor  is  defined  as  the  ratio  of 
the  imaginary  part  to  the  real  part  of  moduli.  The  complex  modulus  is  only 
defined  for  harmonic  excitation.  However,  this  complex  modulus  approach  can  be 
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Figure  2.4  Sketch  of  logarithmic  decrement. 
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used  for  the  analysis  of  a  structural  system  under  arbitrary  excitation  as  any  real 
signal  can  be  represented  as  a  Fourier  series. 

For  a  structure  exhibiting  steady  state  oscillation,  the  loss  factor  r?  can  be 
defined  in  terms  of  energy  quantities  as  follows, 


ri  =  (2.11) 


(InW) 


where  D  is  the  energy  dissipated  per  cycle  and  W  is  defined  as  the  maximum 
energy  stored  during  a  cycle.  Ungar  and  Kerwin  [56]  showed  that  this  definition  of 
W  is  unambiguous  only  in  lightly  damped  systems.  It  was  also  shown  that  for  an 
arbitrary  series-parallel  network  of  viscoelastic  springs,  the  composite  loss  factor 
reduced  to  the  following  simple  expression 

T]  =  ^   (2.12) 

E  K 

where  ry,,  and  W„  are  the  loss  factors  and  the  strain  energy  stored  in  the  n"' 
spring,  respectively. 

Hysteresis  damping  in  a  continuous  composite  system  can  be  modeled  by 
using  the  correspondence  principle  of  viscoelasticity  [57].  According  to  this 
principle,  the  elastic  constants  of  the  materials  are  simply  replaced  by  the 
corresponding  viscoelastic  counterparts.  For  example.  The  dynamic  behavior  of 
such  materials  can  be  defined  in  terms  of  the  complex  modulus 
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E'  =  E'  +  iE"  =  £'(1  +  it])  (2.13) 

where 

E'   =  storage  modulus 
E"  =  loss  modulus 
t]    =  loss  factor 

2.3  Damping  Analysis 
The  analysis  of  a  structure  with  constrained-layer  treatment  is  different 
from  other  situations  because  the  addition  of  the  constraining  and  viscoelastic 
layers  on  the  base  structure  will  change,  in  addition  to  its  mass  and  stiffness 
characteristics,  the  geometry  of  the  structure.  The  stiffness  changes  significantly 
through  the  thickness  of  the  configuration,  varying  from  a  high  value 
corresponding  to  the  base  structure  to  a  low  value  associated  with  the  viscoelastic 
layer.  Ross,  Ungar  and  Kerwin  [58]  developed  what  is  probably  the  most  widely 
used  analysis  to  describe  the  behavior  of  different  type  of  surface  damping 
treatments.  The  fundamental  work  in  the  area  of  the  analysis  of  sandwich 
viscoelastic  beams  was  done  by  Kerwin  [59].  Other  literature  on  this  topic  may  be 
found  in  [60]. 

Since  much  of  the  difficulty  in  designing  constrained-layer  damping 
treatment  is  due  to  complicated  geometries,  it  is  therefore  necessary  to  look  at 
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finite  element  techniques  for  solutions  to  the  problem.  Since  optimization 
techniques  usually  require  a  large  number  of  function  evaluations,  efficient  and 
accurate  finite  elements  are  desirable  to  model  the  composite  configuration. 
Three-dimensional  finite  element  models  of  sandwich  structures  often  result  in 
high-dimensionality  stiffness  matrices  and  associated  significant  computational 
requirements. 

Sun,  Sankar,  and  Rao  [29-33]  have  developed  elements  that  are  specially 
suited  to  modeling  sandwich  structures.  The  finite  element  implementation  of  the 
modal  strain  energy  method  [56,61]  and  the  direct  frequency  response  method 
[61]  are  used  to  analyze  the  free  and  forced  vibration  characteristics  of  the  beam. 
Damping  in  the  system  is  represented  by  using  the  complex  stiffness  approach 
which  derives  from  the  elastic-viscoelastic  correspondence  principle.  In  the  direct 
frequency  response  method,  a  forced  vibration  over  a  range  of  frequencies  is 
considered  and  a  theoretical  response  spectrum  is  generated.  The  loss  factor,  or 
the  overall  damping  of  the  system,  may  then  be  calculated  by  using  the  half- 
power-bandwidth  technique.  The  technique  is  illustrated  in  Figure  2.5.  The  direct 
frequency  response  is  computationally  expensive,  and  it  is  mostly  used  to  verify 
the  predictions  of  other  methods  since  it  needs  no  assumptions  about  the  level  of 
damping.  A  more  frequently  used  method  to  calculate  damping  is  the  modal  strain 
energy  method.  The  method  is  based  on  the  following  eigenvalue  problem, 

[[^  +  [Kg]W)  =  oy'-[Af\{U}  (2.14) 
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Logarithmic  Scale  :  At  3dB  Below  the  Peak  Value  (AT 

RESONANCE  FREQUENCY) 

Linear  Scale  :  At  0.707  the  Peak  Value 
Loss  Factor  =  t\ 
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where 

[K]  =  global  conventional  stiffness  matrix 
[K^]  =  global  geometric  stiffness  matrix 

[M]  =  global  mass  matrix 
iU)    =  global  displacement  vector 

The  natural  frequencies  and  mode  shapes  are  calculated  about  the  initially 
stressed  state  for  the  undamped  system.  The  analysis  can  be  based  on  the 
associated  complex  eigenvalue  problem  or  the  simphfied  real  eigenvalue  problem. 
When  the  analysis  is  performed  based  on  the  simplified  real  eigenvalue  problem, 
the  stiffness  matrix  is  real  and  the  corresponding  nodal  displacements  are  also 
real.  Modal  damping  is  then  calculated  by  this  technique  using  the  modal  strain 
energy  method.  The  loss  factor  is  calculated  as 


It 


rt  =  -i:^   (2.15) 


i-1 

where 

n    =  total  number  of  elements 

r?i  =  loss  factor  of  the  i""  element 

Ej  =  elastic  shear  energy;  equal  to  V2{u}'^[k]{u} 

{u}  =  displacement  vector  of  the  i"'  element 

[k]  =  stiffness  matrbc  of  the  i"'  element 
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This  technique  is  vah'd  only  for  systems  with  relatively  small  levels  of  damping 
where  the  mode  shapes  and  frequencies  of  the  damped  and  the  undamped 
structure  are  similar.  The  analytical  predictions  obtained  by  this  method  were  in 
good  agreement  with  the  experimental  results  of  Mantena  [62]  for  several  damped 
laminated  beam  configurations.  This  analysis  [29-33]  was  used  in  conjunction  with 
optimization  algorithms  for  solutions  of  optimization  problems  associated  with 
viscoelastically  damped  composite  elements  in  the  present  research. 

2.4  Design  Optimization  Formulation 
The  problem  of  optimization  of  a  laminated  composite  beam  for  enhanced 
damping  characteristics  can  be  formulated  as  one  or  all  of  the  following 
subproblems: 

a)  optimal  design  of  the  base  structure. 

b)  design  with  add-on  viscoelastic  damping  treatment. 

c)  integral  design  of  sandwich  beam  consisting  of  base  structure  and 
add-on  damping  treatments. 

Damping  ability  of  a  laminated  composite  beam  can  be  enhanced  by 
optimally  configuring  the  internal  structural  parameters  of  the  composites.  These 
structural  parameters  include  ply  thicknesses,  laminate  thickness,  material 
properties  of  each  laminae,  fiber  orientation,  etc.  By  using  formal  optimization 
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techniques,  the  damping  abihty  can  be  maximized  with  the  satisfaction  of 
prescribed  design  constraints. 

The  use  of  add-on  viscoelastic  damping  is  convenient  and  efficient,  because 
the  configuration  of  the  base  structure  need  not  be  altered.  Given  a  laminated 
composite  beam,  the  treatment  of  viscoelastic  damping  can  be  done  by  attaching 
the  viscoelastic  damping  layers  to  the  surface  of  the  laminated  beam.  Only  one 
viscoelastic  layer  is  required  for  unconstrained-layer  treatments.  One  viscoelastic 
layer  and  a  thin  constraining  layer  are  needed  for  constrained-layer  treatments. 
For  a  two-dimensional  cantilever  beam,  an  arbitrary  length  of  selected  add-on 
damping  tape  can  be  placed  on  the  composite  beam.  Issues  such  as  the  precise 
location  of  damping  material  and  its  quantity  are  resolved  by  expressing  them  as 
decision  variables  in  the  damping-enhancement  optimal  design  problem. 

Placement  of  the  tape  in  a  certain  form  may  create  some  problems  in  the 
analysis.  As  the  tape  moves  from  one  point  to  another  point,  a  new  finite  element 
model  will  be  required,  and  a  new  finite  element  mesh  of  the  structure  is 
necessary.  An  innovative  technique  which  assures  a  smooth  variation  in  node 
locations  during  redesign  is  proposed  in  this  work.  For  a  non-continuous 
placement  of  the  damping  tapes,  a  number  of  unconnected  damping  tapes  can  be 
placed  along  the  length  of  a  beam.  Several  different  ways  of  placing  damping 
tapes  on  a  cantilever  beam  are  shown  in  Figure  2.6.  The  problems  associated  with 
the  discontinuity  of  damping  tapes  manifest  themselves  as  discontinuities  in  the 
node  numbers,  and  the  requirement  of  an  adaptive  preprocessor  which  is  able  to 
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Figure  2.6  Examples  of  damping  treatment  on  cantilever  beam. 
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create  a  new  finite  element  model  for  each  new  structure  created  by  the 
optimizer.  For  a  3-dimensional  laminated  composite  beam,  the  aforementioned 
difficulties  were  overcome  again  with  a  more  sophisticated  adaptive  preprocessor. 
The  mathematical  representation  of  the  presence  or  absence  of  damping 
treatment  on  a  predetermined  section  may  be  resolved  by  assigning  it  a  0/1  type 
of  variable.  In  this  representation,  the  design  problem  contains  discrete  variables, 
and  optimization  algorithms  which  can  account  for  discontinuous  design  variables 
in  an  efficient  manner,  would  be  required. 

The  integral  design  of  the  entire  structure  is  most  effective  for  the 
enhancement  of  the  damping  ability  of  the  structure.  This  method  is  basically  the 
combination  of  the  design  of  the  base  structure  and  the  damping  treatments  at  the 
same  time;  here,  design  constraints  can  be  more  easily  satisfied  due  to  the 
increased  flexibility  in  changing  the  design. 


CHAPTER  3 

NONLINEAR  PROGRAMMING  BASED  OPTIMIZATION  METHODS 


3.1  Introduction 

Most  practical  engineering  optimization  problems  involve  linear  or 
nonlinear  equality/inequality  constraints.  Mathematical  programming  based 
solutions  to  such  problems  can  be  classified  into  three  categories.  The  first 
category  is  that  of  linear  programming  for  which  there  are  relatively  efficient 
solution  algorithms  available.  However,  as  the  name  linear  programming  suggests, 
these  optimization  problems  represent  a  limited  number  of  physical  problems, 
because  the  objective  functions  and  constraints  must  be  strictly  linear.  The  more 
general  optimization  problems  involving  nonlinear  objective  functions  and 
constraints  are  approached  in  two  ways.  The  first  entails  converting  the 
constrained  problems  into  an  unconstrained  optimization  problem,  by  appending  a 
measure  of  constraint  violation  to  the  objective  function.  This  unconstrained 
problem  can  be  solved  by  one  of  the  many  nonlinear  programming  algorithms.  To 
prevent  this  approach  from  becoming  numerically  ill-conditioned,  the  optimal 
solution  is  actually  obtained  by  solving  a  sequence  of  unconstrained  problems, 
corresponding  to  increasing  values  of  the  penalty  on  constraint  violations.  This 
type  of  method  is  referred  to  as  a  sequential  unconstrained  minimization 
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technique  or  SUMT.  The  third  category  consists  of  methods  which  are  commonly 
used  in  practical  engineering  optimization  problems.  This  approach  deals  with 
constraints  directly  in  searching  for  an  optimal  solution,  and  is  classified  as  a 
direct  method.  Such  an  approach  was  adopted  extensively  in  the  present  work, 
primarily  to  obtain  baseline  optimal  solutions  to  test  problems.  These  solutions 
were  then  used  for  comparison  with  solutions  obtained  through  other  schemes 
proposed  in  this  work.  The  direct  method  used  here  is  referred  to  as  the  method 
of  feasible-usable  search  directions,  and  is  summarized  in  the  following  sections. 
A  review  of  mathematical  programming  applications  in  structural  optimization  is 
available  in  [63]. 

3.2  Feasible-Usable  Search  Direction  Method 
An  overview  of  the  optimization  problem  and  its  formal  mathematical 
statement  are  introduced  next.  The  method  of  feasible-usable  search  directions  is 
then  discussed. 

In  an  engineering  optimal  design  problem,  the  quantities  to  be  determined 
are  the  design  variables.  For  a  given  optimal  design  problem,  there  usually  exist 
some  restrictions  dictated  by  the  environment,  processes,  or  resources,  which  have 
to  be  satisfied  in  order  to  produce  an  acceptable  solution.  These  requirements  are 
referred  to  as  constraints,  and  describe  dependencies  among  design  variables  and 
other  parameters.  Constraints  are  generally  expressed  mathematically  by  either 
equalities  or  inequalities.  Among  the  designs  satisfying  all  constraints,  one  has  to 
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be  selected  as  the  best  or  optimal  design  based  on  the  value  of  some  merit  or 
criterion  function.  This  criterion  is  an  integral  part  of  the  optimization  problem 
formulation.  It  is  generally  a  computable  function  of  design  variables  and  other 
parameters,  and  is  called  the  objective  function. 

A  general  constrained  optimization  problem  can  be  mathematically 
formulated  as  follows: 

Minimize:         F{X)  objective  function  (3-1) 
subject  to 

gj{X)  <  0       i=\,m  inequality  constraints  (3.2) 

^'a  (^    "   0       ^=1,^  equality  constraints  (3.3) 

x'  <     <  x"       i=l,p  side  constraints  (3-4) 

where  X  =  [xj,  x,,     Xp]'^  is  a  vector  of  design  variables;  x-  and  x"  represent  the 
lower  and  upper  bounds  of  design  variable  x,,  respectively. 

If  equality  constraints  are  explicit  in  X,  the  number  of  design  variables  can 
often  be  reduced.  Although  side  constraints  (Eqn.  3.4)  can  be  included  in  the 
inequality  constraint  set  (Eqn.  3.2),  it  is  usually  convenient  to  keep  them  separate, 
as  they  define  the  region  of  search  for  the  optimum. 
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Most  mathematical  programming  algorithms  require  an  initial  design  vector 
X^,  where  initial  values  of  each  design  variable  Xj°,  X2°,     Xp°,  are  specified.  From 
this  starting  point,  the  design  is  updated  iteratively  to  an  improved  solution.  The 
most  generally  used  iterative  update  equation  is  as  follow, 

A"?   =  X^'^    +   a^'  5^  (3-5) 

where  superscript  q  is  the  iteration  number,  X  is  the  vector  of  design  variables,  S 
is  the  search  direction,  and  a^'  is  a  scalar  multiplier  determining  the  amount  of 
change  in  X  for  each  iteration.  The  search  direction,  S,  is  first  determined  in  each 
iteration,  and  the  search  moves  in  this  direction  to  update  the  X  vector  according 
to  Equation  3.5. 

The  method  of  feasible  directions,  or  method  of  feasible-usable  search 
directions  [64],  is  introduced  by  an  example  with  two  design  variables  and  two 
inequalities  constraints  with  reference  to  Figure  3.1.  The  figure  shows  the  contour 
lines  of  the  constant  objective  function  F(X),  and  also  shows  the  constraints. 

The  goal  of  this  method  is  to  determine  a  search  direction  which  rapidly 
reduces  the  objective  function  while  maintaining  a  feasible  design.  Consider  a 
design  X°  on  the  constraint  boundary  gi(X)  =  0.  Gradient  vectors  of  the 
objective  function,  VF(X°),  and  the  active  constraint,  Vgi(X°),  are  calculated.  The 
lines  (hyperplanes  in  n  dimensions)  tangent  to  the  contour  of  constant  objective 
and  tangent  to  the  constraint  boundary  are  now  the  linear  approximation  to  the 
problem  at  this  point.  The  next  task  is  to  determine  a  search  direction  S  which 
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Figure  3.1  Illustration  of  the  feasible-usable  search  direction  approach. 
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reduces  the  objective  function  without  violating  the  active  constraint  for  some 
finite  move.  An  S  vector  is  called  a  usable  direction  if  the  objective  function  is 
reduced  or  maintained  at  its  present  value  for  a  finite  move  along  that  direction. 
The  usability  requirement  is  mathematically  defined  as  follows, 

VF{X^)-S  <  0  (3-6) 


A  direction  S  is  called  feasible  if,  for  some  finite  move  in  that  direction, 
the  active  constraint  is  not  violated.  The  feasibility  requirement  can  also  be 
mathematically  defined  as  follows: 

Vg,(Z°)-S<0  (3.7) 


A  region  is  called  feasible-usable  if  any  direction,  S,  inside  that  region 
satisfies  both  Equations  3.7  and  3.8.  The  search  direction  S  is  chosen  to  be 
located  in  the  usable-feasible  region  such  that  the  objective  function,  for  a  small 
move  step,  is  reduced  while  the  feasibility  is  maintained. 

Since  a  method  such  as  the  feasible-usable  direction  approach  is  essentially 
a  local  search  technique,  only  the  relative  optimum  will  be  attained  if  the  design 
space  is  nonconvex.  A  design  can  be  treated  as  the  globally  optimal  design  with  a 
higher  degree  of  confidence  only  if  a  large  number  of  optimization  processes  with 
different  starting  values  of  vector  X  converge  to  the  same  final  design.  If  a 
selected  number  of  optimization  processes  results  in  different  "optimal"  designs,  it 
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usually  indicates  that  the  design  space  is  multimodal  or  nonconvex.  This  problem 
can  generally  be  solved  by  either  switching  to  more  sophisticated  algorithms  or  by 
considering  more  initial  designs  to  increase  the  probability  of  locating  the  global 
optimal  design.  It  should  be  noted  that  any  optimization  method,  except  the 
exhaustive  search  technique,  provides  no  guarantee  of  locating  the  global 
optimum. 

To  improve  computational  efficiency  in  nonlinear  optimization  problems,  a 
piecewise  linearization  of  the  design  space  is  frequently  adopted.  In  this  approach, 
the  nonlinear  programming  problem  (both  objective  function  and  constraints)  is 
first  linearized  about  the  starting  design.  The  solution  to  this  approximate  linear 
problem  is  obtained  within  some  prescribed  move  limits  so  as  to  maintain  the 
integrity  after  linear  approximation.  The  problem  is  next  linearized  about  the  new 
solution,  and  the  process  repeated  till  convergence. 

3.3  Modified  Branch-and-Bound  Search 
The  discrete  nature  of  the  design  variables  in  most  structural  optimization 
problems  has  only  received  limited  attention.  Solutions  to  such  problems  are 
relatively  more  difficult  to  obtain  than  for  optimization  problems  with  only 
continuous  design  variables.  Discrete  design  variables  in  structural  optimization 
have  often  been  treated  as  continuous  design  variables.  Upon  obtaining  the 
optimal  design  for  this  simplified  problem,  design  variables  which  represented 
discrete  design  variables  were  then  adjusted  to  the  nearest  discrete  value.  This 
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simplified  approach  worked  well  for  optimization  problems  with  discrete  design 
variables  whose  values  were  spaced  reasonably  close  to  one  another.  However, 
this  simple  rounding  procedure  often  failed  in  locating  the  true  optimal  solution, 
or  resulted  in  infeasible  solutions  [65].  One  branch  of  mathematical 
programming,  referred  to  as  integer  programming,  offers  a  formal  solution  for 
optimization  problems  using  integer/discrete  design  variables.  Early  work  in 
obtaining  systematic  solutions  to  the  integer  linear  programming  problem  are 
described  by  Gomory  [66].  The  branch-and-bound  algorithms  that  emerged  later 
are  based  on  the  enumeration  of  the  space  of  the  feasible  integer  solutions.  A 
more  detailed  explanation  of  the  branch-and-bound  approach  is  available  in  [67]. 

The  general  framework  for  solving  an  integer  programming  problem 
involves  decomposing  the  original  problem  into  subproblems,  modifying 
constraints  to  enlarge  feasible  domains,  and  finally  a  process  referred  to  as 
fathoming.  Fathoming  involves  checking  a  solution  for  feasibility,  and  establishing 
optimality. 

A  modified  mixed  integer  and  discrete  programming  algorithm  was 
proposed  by  Hajela  and  Shih  [68].  The  algorithm  employs  a  strategy  where  a 
systematic  search  of  continuous  solutions  is  made,  in  which  the  discrete  and 
integer  variables  are  successively  forced  to  assume  specific  values.  The  basic 
solution  strategy  used  in  this  method  for  nonlinear  mixed  integer  programming 
problem  is  a  variant  of  the  approach  proposed  by  Garfinkel  and  Nemhauser  [69]. 
The  logical  structure  of  the  set  of  solutions  was  constructed  as  a  binary  tree.  A 
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modified  feasible  directions  algorithm  [70]  was  used  in  the  solution  of  the 
continuous  nonlinear  programming  problem,  with  piecewise  linear  representation 
of  the  objective  function  and  constraints. 

3.4  Min-Max  Global  Criterion  Approach 
Decisions  in  engineering  design  typically  require  allocation  of  resources  to 
satisfy  multiple,  and  frequently  conflicting,  requirements.  For  instance,  in  the 
design  of  a  viscoelastically  damped  cantilever  beam,  minimum  weight  and 
maximum  damping  are  two  criteria  which  always  conflict  with  each  other. 
Nevertheless,  these  two  criteria  are  both  important  characteristics.  Adding  on 
viscoelastic  damping  tape  would  increase  the  damping  ability  of  the  beam  while, 
at  the  same  time,  increase  the  weight.  One  approach  to  solving  this  multicriterion 
optimization  problem  is  to  select  weight  as  the  objective  function,  and  to  use 
damping  as  a  design  constraint,  or  vice  versa.  The  apparent  simplicity  afforded  by 
this  method  is  attractive;  however,  it  should  be  noted  that  the  treatment  of 
criteria  as  constraints  does  not  yield  the  same  optimal  design  as  would  be 
obtained  when  solving  the  optimum  design  problem  as  one  possessing  multiple 
objectives  [28].  In  the  above  optimization  example,  the  influence  of  the  relative 
importance  of  structural  weight  or  modal  damping  on  the  distribution  of 
viscoelastic  damping  material  can  be  obtained  only  by  formally  solving  a  two- 
criterion  problem  with  different  weighting  coefficient  combinations  on  criteria. 
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In  order  to  help  understand  the  solutions  to  a  multicriterion  optimization 
problem,  a  frequently  used  expression,  "Pareto  optimal  solution",  needs  to  be 
defined.  Given  two  designs  Xi(xi,X2,...x„)  and  Xj(xi,X2,...Xn,)  in  an  n-criterion 
minimization  problem,  the  design  Xj  is  referred  to  as  being  partially  less  than  the 
design  X^  if  each  element  of  the  criterion  vector  of  the  design  Xj  is  less-than-or- 
equal  to  the  corresponding  element  of  design  Xj,  and,  for  at  least  one  element, 
the  relation  is  strictly  less-than.  If  X|  is  partially  less  than  X^,  X,  is  said  to 
dominate  X^,  or  Xj  is  inferior  to  Xj.  In  a  given  set  of  designs,  a  population  of 
designs  which  are  not  dominated  by  any  other  design  in  the  given  set  is  said  to  be 
non-dominated  or  non-inferior.  The  set  of  non-dominated  or  non-inferior  designs 
are  referred  to  as  the  Pareto-optimal  designs  which  are  referred  to  as  Pareto- 
optimal  solutions  in  an  optimization  problem. 

An  effective  min-max  variant  of  the  global  criterion  approach  for 
multicriterion  optimization  problems  was  introduced  by  Hajela  and  Shih  [71].  It  is 
referred  to  as  a  global  criterion  method  [72],  and  is  broadly  classified  as  belonging 
to  a  category  of  solution  methods  with  no  articulation  of  preference.  In  such  an 
approach,  a  metric  function  is  formulated  to  represent  the  distance  between  the 
ideal  solution  and  the  optimum  solution,  and  a  minimization  of  this  function 
results  in  the  true  optimum.  The  ideal  solution  is  an  n-dimensional  vector  for  an 
optimization  problem  with  n  design  criteria.  An  optimal  value  for  each  individual 
criterion  defines  a  component  of  the  n-dimensional  vector.  These  'n'  components 
of  the  ideal  solution  vector  are  obtained  by  considering  the  optimization  problem 
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for  each  criterion  separately.  An  application  of  the  global  criterion  method  is 
illustrated  by  Figure  3.2.  This  figure  shows  the  feasible  and  infeasible  sectors  of 
the  space  of  two  criterion  functions.  If  the  ideaJ  solution,  f^(X),  was  located  inside 
the  feasible  region,  no  further  effort  would  be  necessary.  However,  the  ideal 
solution  is  typically  infeasible,  and  the  task  becomes  one  of  finding  the  feasible 
solution  which  is  closest  to  the  ideal  solution.  The  Pareto  optimal  solution  is 
obtained  in  such  a  fashion  that  any  change  of  the  optimal  design  vector  will  cause 
adverse  effects  on  at  least  one  of  criteria. 

One  can  define  a  vector  objective  function  f(X)  dependent  on  the  design 
variable  vector  X  in  the  following  fashion: 

m  =  [f,(X}/,{X),...f^{x),.../,(x)f 

(3.8) 

where  f,(X)  is  the  i"'  objective  function  or  criterion.  Likewise,  {^(X"^)  is  the  ideal  or 
optimal  value  if  the  i""  objective  function  is  considered  as  the  sole  criterion  of  a 
single-criterion  optimization  problem,  and  X''^  is  the  optimal  design  vector  for 
such  a  criterion. 

Minimization  of  the  metric  distance  between  the  infeasible  ideal  solution 
and  the  candidate  feasible  solution  by  a  min-max  variant  of  the  global  criterion 
approach  leads  to  the  optimization  problem  described  by  the  following 
mathematical  statement. 
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Figure  3.2  Graphical  representation  of  Pareto  (noninferior)  solutions  in  two-criterion 
minimization  problem. 
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Min  Max 


(3.9) 


Solution  of  the  above  optimization  problem  yields  the  best  compromise 
solution  in  which  all  criteria  are  considered  equally  important.  Use  of  weighting 
coefficients  can  be  introduced  in  conjunction  with  this  method  to  rank  the 
importance  of  each  criterion.  The  min-max  problem  can  be  restated  as  follows: 


where  W;  is  the  weighting  coefficient  representing  the  relative  importance  of  the  i"" 
criterion. 

A  scalar  variable  y3  is  introduced  to  transform  the  min-max  problem  of 
(3.10)  into  an  equivalent  scalar  optimization  problem.  This  single-criterion 
optimization  problem  is  as  follows: 


Min  Max  w- 


i=l,2,...,k 


(3.10) 


Minimize 


(3.11) 


subject  to  the  following  additional  constraints: 
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f,iX)-f^{X") 


ax-') 


-  p  <  0 


i=\,2,-,k 


(3.12) 


In  a  situation  where  the  design  variable  set  is  a  mixture  of  integer,  discrete, 
and  continuous  design  variables,  it  can  be  shown  [71]  that  the  ideal  solution 
should  be  selected  as  the  one  obtained  by  treating  all  design  variables  as 
continuous. 


3.5  Numerical  Implementation 
The  modified  branch-and-bound  method,  in  conjunction  with  the  min-max 
variant  of  the  global  criterion  approach,  was  implemented  in  determining  the 
multicriterion  optimal  design  of  a  sandwich  cantilever  beam  to  enhance  its 
damping  ability.  Damping  ability  in  a  structural  composite  or  metal  beam  can  be 
enhanced  passively,  via  an  external  viscoelastic  damping  layer  treatment.  The 
design  space  for  this  problem  consists  of  both  continuous  and  discrete/integer 
design  variables.  An  aluminum  cantilever  beam  serves  as  a  base  structure  of  the 
structural  system  for  which  maximum  damping  in  several  modes  is  desired. 
Constrained  viscoelastic  layers  were  added  to  increase  the  damping  ability  of  the 
structure.  The  structural  system  and  two  methods  of  adding  damping  layers  are 
shown  in  Figures  3.3  and  3.4.  Material  properties  for  this  structural  system  are 
given  in  Table  3.1.  In  the  first  case,  the  add-on  viscoelastic  tape  was  applied 
starting  from  the  root  of  the  cantilever  beam  with  the  total  length  of  the  tape 
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Figure  3.3  Sandwich  beam  structure. 
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Figure  3.4  Sandwich  beam  structure. 
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Table  3.1  Material  properties  of  the  structural  system. 


cuiisudining  layer,  ano 
base  structure 

viscoeiastic  layer 

Young's  modulus 
(GN/m^) 

69.0 

0.0021 

Poisson's  ratio 

0.3 

0.499 

Damping  factor 

0.0 

0.10 

Density 
(kg/m^) 

2700.0 

970.0 
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indicated  by  X5  (X4  =  0).  In  the  other  case,  the  add-on  damping  tape  was  applied 
starting  from  an  arbitrary  point  on  the  beam.  The  distance  between  the  root  of 
the  beam  and  the  starting  point  is  indicated  by  X4,  and  the  total  length  of  the  tape 
is  indicated  by  X5. 

One  integrated  approach  to  the  optimal  design  of  the  viscoelastically 
damped  structure  is  to  configure  both  the  basic  structure  and  add-on  damping 
layers.  The  optimal  configuration  is  then  determined  by  the  thickness  of  the 
viscoelastic  layer,  the  thickness  of  the  constraining  layer,  and  the  thickness  of  the 
base  structure.  They  are  denoted  as  Xj,  X2,  and  X3,  respectively.  The  remaining 
two  design  variables,  X^  and  X5,  result  from  how  the  damping  layers  were  applied. 

Finite  element  analysis  was  used  to  evaluate  damping  and  natural 
frequencies,  as  well  as  the  static  displacement  in  the  structure  resulting  from  the 
different  applied  lengths  of  the  constrained  viscoelastic  layer.  The  structure  of  the 
finite  element  model  representing  the  three  layer  sandwich  is  shown  in  Figure  3.5. 

The  base  structure  and  constraining  layer  were  modeled  by  using  a 
specially  developed  three-node,  seven-degree-of-freedom,  offset  beam  element  as 
shown  in  Figure  3.6.  The  element  contains  shear-deformable  characteristics,  which 
are  important  in  modeling  fiber-reinforced  composites.  A  key  feature  of  this 
model  was  its  accountability  of  the  coupling  between  the  stretching  and  bending 
deformations.  This  allowed  the  nodes  to  be  offset  to  one  surface  of  the  beam, 
coincident  with  the  nodes  of  the  adjoining  element.  The  viscoelastic  core  was 
modeled  using  a  rectangular  plane  stress  element  that  is  compatible  with  the 
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Figure  3.5  Typical  finite  element  model  of  the  sandwich  beam. 
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Figure  3.6  Offset  beam  element. 
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offset  beam  element.  The  stiffness  of  the  offset  element  was  obtained  from  a 
strain  energy  formulation  which  was  based  on  the  following  assumed  displacement 
field: 


u{x^)=u,^(z-!L)y\f{x)  (3.13) 


o)(x^)  =  w{x)  (3-14) 
i^{x^)  =  ^^{x)  (3.15) 


where  Uq  and  i|r  were  obtained  by  linear  interpolation  and  a  quadratic 
interpolation  was  used  to  compute  o).  A  consistent  mass  matrix,  evaluated  from 
the  kinetic  energy,  was  used  in  this  analysis.  A  detailed  description  of  the  offset 
beam  element  is  given  in  [31].  The  loss  factor  of  the  system  was  determined  by 
the  strain  energy  based  approach  discussed  in  the  previous  chapter. 

The  vector  of  design  objectives  included  minimum  structural  weight  and 
maximum  loss  factors  in  the  first  three  modes  of  vibration.  The  constraint  set 
comprised  lower  and  upper  bounds  on  the  first  three  natural  frequencies,  and  the 
maximum  displacement  when  a  static  load  of  1-N  was  applied  at  the  tip  of  the 
beam.  These  are  summarized  in  Table  3.2.  Design  variables  Xj  through  X5  are 
shown  in  Figures  3.3  and  3.4.  Thicknesses  of  the  constraining  and  viscoelastic 
layers,      and  X,,  were  considered  discrete  variables,  with  admissible  values 
selected  from  the  set  shown  in  Table  3.2.  All  other  design  variables  were  assumed 
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Table  3.2  Design  constraints  and  allowable  discrete  values  of  design  variables  Xi 
and  X,. 


Xj,  X2 

0.0001  0.0002  0.0003  0.0004  0.0005  0.0006 
0.0007  0.0008  0.0009  0.0010  0.0010  0.0012 
0.0013  0.0014  0.0015  0.0016  0.0017  0.0018 
0.0019  0.0020  (m) 

Design 
Constraints 

20.0  <  freqi  <  50.0  (Hz) 
130.0  <  freq,  <  200.0  (Hz) 
400.0  <  freqj  <  700.0  (Hz) 
tip  displacement  <  0.001  m 
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continuous.  The  mathematical  statement  of  the  optimum  design  problem  can  be 
written  as  follows: 

Minimize  (-F^,-F,,-F„F^) 

(3.16) 

Subject  to       g.  <  0,  j=l,2,...m 

Here  Fj,  F,,  and  F3  represent  the  loss  factors  of  the  first  three  vibration  modes,  F^ 
is  the  structural  weight,  and  g^ 

are  the  m  inequality  constraints  for  the  problem  as  specified  in  Table  3.2.  Two 
distinct  cases  were  considered  by  setting  X4  to  zero,  and  by  allowing  X4  to  assume 
finite  values,  respectively. 

The  locations  and  lengths  of  the  damping  tape  were  adjusted  using  an 
adaptive  preprocessor.  Ten  elements  along  the  length  of  the  base  beam  were  used 
in  the  finite  element  analysis  throughout  the  optimization  process.  A  uniform 
node  spacing  was  used  in  the  situation  that  X4  and  X5  are  both  multiples  of  one- 
tenth  of  the  length  of  the  beam.  When  the  length  of  an  untreated  free  or  fixed 
end  is  small,  an  unsmooth  node  spacing  may  result  when  an  equal-length-element 
representation  of  the  treated  length  is  used.  To  avoid  this  situation,  stretching 
techniques  were  utilized  to  assure  a  smooth  variation  in  node  spacing.  A 
minimum  node  spacing  ratio  for  any  two  consecutive  finite  elements  so  obtained 
resulted  in  an  improved  level  of  accuracy  on  the  numerical  prediction.  Figure  3.7 
shows  a  beam  illustrating,  from  top  to  bottom,  uniform  node  spacing,  unstretched 
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Figure  3.7  Illustrations  of  node  spacing  arrangements. 
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node  spacing,  and  stretched  spacing.  The  adaptive  preprocessor  also  possesses 
automatic  node  numbering  ability  to  automatically  create  a  finite  element  model 
based  on  the  given  set  of  design  variables. 

The  min-max  global  criterion  approach  was  used  to  transform  the 
optimization  statement  (3.16)  as  follows: 


Minimize 


(3.17) 


subject  to  m  original  constraints  plus  the  following  four  additional  constraints: 


w, 


f,{X)-f,{X'') 


f,iX)-f,{X") 


f,{X)-f^{X'') 


m-ux^') 


13  <  0 


-  13  <  0 


-  (3  <  0 


13  <  0 


(3.18) 


where  W^,  i  =  l,4,  are  weighting  parameters  and  each  weight  is  used  to  emphasize 
the  relative  importance  of  the  i"^  criterion.  The  nonlinear  optimization  problem 
presented  in  (3.17  and  3.18)  was  solved  by  the  feasible-direction  method  with  the 
modified  branch-and-bound  technique  [68]. 
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For  case  1,  the  separately  attainable  optimal  designs  of  each  component  of 
the  objective  function  are  shown  in  Table  3.3.  The  design  variables  were  all 
assumed  continuous  in  finding  the  ideal  solutions.  Assigning  different  weighting 
coefficients  to  the  components  of  the  objective  function  permitted  the  generation 
of  a  series  of  non-inferior  optimal  solutions,  as  presented  in  Table  3.4.  These  non- 
inferior  solutions  were  such  that  no  component  of  the  objective  vector  could  be 
improved  without  adversely  affecting  another  component.  Similar  solutions  for 
case  2  are  presented  in  Tables  3.5  and  3.6.  The  inclusion  of  the  X4  variable  results 
in  no  improvement  in  F,'",  Fj'",  and  F^^  However,  F,*''  (loss  factor  in  second 
mode)  increased  by  2%  and  the  length  of  the  damping  layer  X5  decreased  by 
18%.  Permitting      to  assume  finite  values  yields  a  noticeable  improvement  in  the 
objective  function,  only  when  the  weighting  coefficients  associated  with  F3 
outweigh  those  of  Fj  and  F3.  This  was  because  the  middle  region  of  the  beam 
undergoes  the  most  shear  deformation  for  a  displacement  corresponding  to  the 
second  mode  of  vibration.  With  X4  allowed  to  vary,  and  for  weighting  coefficients 
set  as  (0.0,  0.5,  0.0,  0.5),  the  damping  factor  of  the  second  mode  was  increased  by 
5%,  while  the  weight  of  the  structure  was  reduced  by  1%. 
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Table  3.3  Ideal  solutions  for  case  1. 


Fi"*^^'  =  0.01825 

F,""^"  =  0.02003 

Fj"'^^'  =  0.01833 

F;'"=^'  =  1.5860 

X"'"'(Fj) 

-  (  0.0020,  0.0001,  0.00454,  0.262  ) 

X"^"'(F2) 

=  (  0.0020,  0.0001,  0.00446,  0.358  ) 

=  (  0.0020,  0.0001,  0.00445,  0.400  ) 

=  (  0.0001,  0.0001,  0.00499,  0.000  ) 
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Table  3.4  Optimal  solutions  for  problem  in  case  1. 


\A/ j^i  rthti'nrT 

vv  cigjiiing 

uesign  vanaDies 

Loss  Factor 

weignt 

^^omomdiion 

V    V    V  V 

Mode  1 

W  J- W4 

(m) 

Mode  2 

MnHp  3 

70 

n  nni/i 

n  MAO 

U.UUUl 

U.Ui!)4 

.zu 

U.UU44 

0.013  / 

.40 

0.359 

.40 

0.0019 

0.0172 

1.9350 

.30 

0.0001 

0.0174 

.10 

0.0046 

0.0098 

.20 

0.309 

.20 

0.0018 

0.0167 

1.9198 

.40 

0.0001 

0.0187 

.20 

0.0044 

0.0156 

.20 

0.357 

.00 

0.0015 

0.0149 

1.8397 

.25 

0.0001 

0.0163 

.25 

0.0044 

0.0140 

.50 

0.357 

.00 

0.0015 

0.0149 

1.8382 

.50 

0.0001 

0.0163 

.00 

0.0044 

0.0139 

.50 

0.355 
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Table  3.5  Ideal  solutions  for  case  2. 


Fj"^"'  =  0.01825 


=  0.02037 


F/^"'  =  0.01833 


F,"""-'^  =  1.5829 


X.deai^p^^  =  (  0.0020,  0.0001,  0.00454,  0.000,  0.262  ) 


X"''=^'(F2)  =  (  0.0020,  0.0001,  0.00446,  0.060,  0.294  ) 


X"'"'(F3)  =  (  0.0020,  0.0001,  0.00445,  0.000,  0.400  ) 


X""^^'(F,)  =  (  0.0001,  0.0001,  0.00498,  0.054,  0.000  ) 
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Table  3.6  Optimal  solutions  for  problem  in  case  2. 


Weighting 

Design  Variables 

Loss  Factor 

Weight 

Combmation 

Xj,X-,,X3,X4,X5 

Mode  1 

W1-W4 

(m) 

Mode  2 

Mode  3 

.20 

0.0017 

0.0115 

1.8765 

.20 

0.0001 

0.0172 

.20 

0.0045 

0.0129 

.40 

0.022 

0.319 

.40 

0.0019 

0.0172 

1.9350 

.30 

0.0001 

0.0174 

1  A 

.10 

0.0046 

0.0098 

.20 

0.000 

0.309 

.20 

0.0018 

0.0167 

1.9198 

.40 

0.0001 

0.0187 

.20 

0.0044 

0.0156 

.20 

0.000 

0.357 

.00 

0.0015 

0.0083 

1.8433 

.25 

0.0001 

0.0155 

.25 

0.0045 

0.0144 

.50 

0.038 

0.333 

.00 

0.0017 

0.0066 

1.8288 

.50 

0.0001 

0.0171 

.00 

0.0046 

0.0098 

.50 

0.059 

0.262 

CHAPTER  4 
GENETIC  SEARCH  METHODS 


4.1  Introduction 

Mathematical  nonlinear  programming  algorithms  have  emerged  as  the 
method  of  choice  for  applications  in  engineering  optimization  problems.  They 
provide  a  general  approach  for  obtaining  solutions  to  both  single  and 
multi-objective  design  problems  with  a  mix  of  equality  and  inequality  constraints. 
The  more  efficient  of  this  class  of  methods  are  generally  gradient  based,  and 
require  at  least  the  first-order  derivatives  of  both  objective  and  constraint 
functions  with  respect  to  the  design  variables.  With  this  "slope-tracking"  ability, 
gradient-based  methods  can  easily  identify  a  relative  optimum  closest  to  the  initial 
guess  of  the  optimum  design.  There  is  no  guarantee  of  locating  the  global 
optimum  if  the  design  space  is  known  to  be  nonconvex  [73].  These  methods  are 
also  inadequate  in  problems  where  the  design  space  is  discontinuous,  because  the 
derivatives  of  both  the  objective  function  and  constraints  may  become  singular 
across  the  boundary  of  discontinuity. 

In  engineering  design  problems,  the  mix  of  continuous,  discrete,  and 
integer  design  variables  has  been  approached  by  treating  all  variables  as 
continuous,  and  then  rounding  specific  variables  either  up  or  down  to  the  nearest 
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integer  or  discrete  variable.  Branch-and-bound  techniques  based  linear 
programming  or  combinatorial  programming  offers  a  formal  solution  to  this  class 
of  optimization  problems  with  a  mix  of  integer,  discrete  and  continuous  design 
variables.  This  strategy  consisted  of  a  systematic  search  of  continuous  solutions  in 
which  the  discrete  and  integer  variables  were  successively  forced  to  assume 
specific  values.  However,  in  doing  so,  the  original  optimization  problem  was 
undesirably  expanded  to  a  large  number  of  sub-optimization  problems. 

Exhaustive  search  and  random  search  methods  are  among  the  simplest  and 
most  robust  strategies  for  automated  optimum  design  problems.  These  methods 
can  work  on  almost  all  kinds  of  design  spaces  and  without  any  restriction  on  types 
of  design  variables.  An  improvement  on  the  simple  enumerative  techniques  is 
possible  with  methods  such  as  random  walk  and  random  walk  with  direction 
exploitation.  The  only  drawback  is  that  these  methods  often  require  thousands  of 
function  evaluations  to  achieve  the  optimum,  even  for  the  simplest  of  problems.  It 
is  hence  crucial  to  examine  alternative  strategies  for  optimal  structural  design 
problems,  which  need  less  computational  effort  than  required  by  the  enumerative 
search  techniques,  and  are  also  not  susceptible  to  convergence  to  a  local  optimum 
as  exhibited  by  gradient-based  nonlinear  programming  algorithms.  Genetic 
algorithms  (GAs)  as  proposed  by  Holland  [34],  have  the  potential  to  successfully 
fill  this  gap. 

Genetic  algorithms  belong  to  a  category  of  stochastic  search  techniques, 
where  only  the  most  promising  regions  of  the  design  space  are  enumerated  to 
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locate  the  optimal  design.  These  algorithms  have  their  philosophical  basis  in 
Darwin's  theory  of  survival  of  the  fittest.  Analogous  to  the  natural  process  where 
a  population  of  a  given  species  adapts  to  a  natural  habitat,  a  population  of  designs 
is  created  and  is  then  allowed  to  adapt  to  the  design  requirements.  Designs  that 
do  not  adapt  in  a  favorable  manner  to  the  requirements  are  eliminated  from 
consideration.  The  mechanism  of  adaptation  borrows  extensively  from  principles 
of  biological  evolution,  in  that  basic  characteristics  of  designs  in  one  population 
are  transferred  to  a  population  in  another  generation  through  gene  transfer 
operators.  Stated  differently,  design  alternatives  representing  a  population  in  a 
given  generation  are  allowed  to  reproduce  and  cross  among  themselves  with  bias 
allocated  to  the  most  fit  members  of  the  population.  Combination  of  the  most 
favorable  characteristics  of  the  mating  members  of  the  population  results  in  a 
progeny  population  that  is  more  fit  than  the  parent  population.  If  the  measure 
which  indicates  the  fitness  of  a  generation  is  also  the  desired  goal  of  a  design 
process,  successive  generations  produce  better  values  of  the  objective  function. 

The  mechanics  of  genetic  search,  though  simple  to  implement,  encompass 
features  that  render  the  approach  highly  applicable  to  the  problem  of  search  in  a 
nonconvex/disjoint  design  space  with  a  mix  of  continuous,  discrete,  and  integer 
design  variables.  These  desirable  characteristics  are  largely  attributed  to  the  fact 
that  genetic  search  moves  from  a  population  of  designs  to  another  population  of 
designs;  this  is  in  contrast  to  the  point  to  point  search  available  in  traditional 
mathematical  programming  methods,  and  therefore  offers  a  better  possibility  of 
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locating  a  global  optimum.  Furthermore,  genetic  algorithms  work  on  a  coding  of 
the  design  variables  rather  than  the  variables  themselves.  This  allows  for  an 
efficient  treatment  of  integer  and  discrete  variables.  The  terminology  of  genetic 
search  and  its  principal  components  are  summarized  next. 


4.2  Genetic  Algorithms 
The  basic  approach  in  genetic  algorithms  is  to  represent  possible  solutions 
to  a  given  problem  by  a  population  of  bit  strings  of  finite  lengths,  and  to 
subsequently  use  transformations  analogous  to  biological  reproduction  and 
evolution  to  improve  and  vary  the  coded  solutions.  In  natural  populations,  genetic 
information  stored  in  chromosomal  strings  evolves  over  generations  to  adapt 
favorably  to  a  static  or  changing  environment.  The  chromosomal  structure 
represents  a  generational  memory,  and  is  altered  through  chromosomal  string 
inversions,  occasional  mutation,  and  a  crossover  of  genetic  information  between 
reproducing  members.  In  an  elitist  reproduction  strategy,  those  members  of  the 
population  that  are  deemed  most  fit  are  selected  for  reproduction,  and  are  given 
the  opportunity  to  strengthen  the  chromosomal  makeup  of  the  progeny 
generation.  This  approach  is  facilitated  by  defining  a  fitness  function  or  a 
measure  indicating  the  "goodness"  of  a  member  of  the  population  in  a  given 
generation  during  the  evolution  process.  For  unconstrained  maximization 
problems,  the  objective  function  could  serve  as  the  fitness  function.  The  inverse  of 
the  objective  function,  or  the  difference  between  a  large  number  (say  the 
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maximum  value  of  the  objective  function  in  the  population)  and  the  objective 
function  value  for  each  member  can  be  used  as  the  fitness  measure  in  a 
minimization  problem.  For  constrained  optimal  design  problems,  an  exterior 
penalty  function  formulation  can  be  adapted  to  transform  a  constrained 
optimization  problem  into  an  unconstrained  one. 

In  biological  systems,  generational  memory  is  preserved  and  transferred  to 
progeny  populations  in  the  form  of  chromosomes.  A  number  of  chromosome 
strings  comprises  a  genotype,  the  total  genetic  makeup  for  an  organism.  In  order 
to  use  genetic  algorithms,  an  artificial  chromosome-like  string  must  be  devised  to 
represent  the  n-dimensional  solution  space.  One  simple  yet  effective  approach  to 
accomplish  this  is  to  represent  each  design  variable  by  a  finite  length  binary  string 
and  then  connect,  head-to-tail,  all  n  strings  into  a  single  binar>'  string.  Various 
simulations  of  genetic  evolution  and  adaptation  are  conceivable.  Three  principal 
components  of  the  gene-transformation  mechanism  in  this  artificial  evolution  and 
adaptation  simulation  are  reproduction,  crossover,  and  mutation.  These 
transformations  are  best  described  with  reference  to  a  specific  optimization 
problem. 

Consider  the  following  function  minimization  problem: 


Minimize 


F{x^^^)  =  xf  +X^X2  +X2 


(4.1) 


Subject  to:         =  x^+xfx^  <  0 


(4.2) 


(4.3) 
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xr  ^  X,  <  xr  i  =  1,2  (4.4) 


To  use  genetic  search,  the  constrained  minimization  problem  is  first  converted 
into  an  unconstrained  problem  using  the  exterior  penalty  function  formulation, 
resulting  in  the  following  problem: 

Minimize  F'  =  F  +  p  (4-5) 

where  P  is  the  penalty  parameter,  and  is  typically  given  by, 

P=rj:<g^>^  (4.6) 

where  r  is  a  penalty  parameter  of  the  form  encountered  in  the  exterior  penalty 
function  approach  [74],  and  <gj>  represents  the  violated  constraints. 

Careful  consideration  must  be  given  to  the  selection  of  the  penalty  function 
P.  Numerical  experiments  have  shown  that  genetic  search  with  a  penalty  function 
formulated  as  in  Eqn.  (4.6),  exhibits  slow  convergence.  This  was  largely  due  to  the 
fact  that  a  severely  violated  constraint  yields  a  value  of  P  which  overwhelms  the 
objective  function  value  F.  A  more  detailed  explanation  of  this  slow  convergence 
will  be  given  in  a  later  section.  In  order  to  overcome  this  situation,  the  following 
bounding  strategy  was  adopted.  If  the  average  fitness  of  feasible  designs  is  F. . 

O  civ ' 

then  a  limiter  value  of  the  penalty  L  is  selected  as. 
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L  =k*F„ 


(4.7) 


where  k  is  of  order  2;  the  penalty  that  is  appended  to  an  infeasible  design  is  then 
obtained  as  follows: 


otherwise 


(4.8) 


The  effect  of  this  scaling  is  to  prevent  radical  departures  in  the  value  of 
the  penalty  term  from  the  specified  value.  If  the  slope  parameter,  a,  is  0.0, 
penalty  for  all  violated  designs  with  P  >  L  is  set  to  L  If  instead,  a  is  assigned  a 
small  value  of  the  order  0.1,  then  the  extent  of  constraint  violation  due  to  severely 
violated  and  less  violated  designs  varies  linearly  from  L,  albeit  with  a  small  slope. 
To  convert  this  function  minimization  into  a  fitness  maximization  as  required  by 
GA's,  the  following  fitness  function  is  created. 


f.  =  F*    -  F* 

J I         max       ^  ( 


(4.9) 


where  f,  is  the  fitness  of  the  i"^  design,  and  F'„,,^  is  the  maximum  value  of  F*.  The 
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F      can  be  chosen  as  the  maximum  F'  in  the  population  of  that  given 
generation,  or  of  several  previous  generations.  A  fixed  large  value  of  F\^^^  can  be 
used  throughout  the  entire  genetic  search  provided  that  there  is  not  any  F*  which 
is  larger  than  F'^^^^  However,  genetic  search  with  a  variable  value  of  F'^^^. 
selected  to  be  slightly  larger  than  the  largest  F'  in  the  present  or  including  past 
few  generations,  usually  converges  faster  than  with  a  fixed  F'„,3^.. 

To  obtain  a  bit  string  representation  of  the  design  variables  in  this 
problem,  each  x,  can  be  converted  into  a  binary  string  of  O's  and  I's.  For 
purposes  of  discussion,  we  choose  a  5  digit  binary  number,  with  the  minimum  and 
maximum  values  of  x^  denoted  as  x,"^'"  and  x,"''^  respectively,  represented  by  the 
following  binary  numbers. 


x,"""  =  00000 
x,"'"'  -  11111 


A  linear  scaling  can  be  introduced  to  convert  intermediate  values  of  the  binary 
number  into  design  variable  values.  The  binary  string  representations  for  Xj  and 
X,  can  be  placed  head-to-tail  to  create  a  10-digit  number,  also  referred  to  as  a 
schema,  which  represents  a  solution  to  the  problem.  Several  such  10-digit  binary 
strings  are  defined  to  constitute  a  population  of  designs,  which  includes  a  mix  of 
feasible  and  infeasible  designs.  The  fitness  f,  corresponding  to  each  member  of  the 
population  is  computed  before  invoking  the  genetic  transformations. 
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4.2.1  Reproduction 

The  simulation  of  genetic  evolution  here  is  contrived  in  that  the  population 
size  is  forced  to  remain  unchanged  during  the  evolution  process;  two  mating 
parents  create  only  two  progenies  but  are  themselves  eliminated.  The 
reproduction  process  is  one  that  biases  the  search  toward  producing  more  fit 
members  in  the  population  and  ehminating  the  less  fit  ones.  One  simplistic 
approach  to  selecting  members  from  an  initial  population  to  participate  in  the 
reproduction  is  to  assign  each  member  a  probability  of  selection  on  the  basis  of  its 
fitness.  If  fj  is  the  fitness  measure  of  the  i*"  member,  the  probabiHty  for  the  i"^ 
member  to  participate  in  reproduction  for  each  of  m  selections  can  be  assigned 
as: 


m 


(4.10) 


where  m  is  the  population  size.  After  m  selections,  a  new  population  pool  is 
created.  The  new  pool  is  the  same  size  as  the  original  pool,  but  has  a  higher 
average  fitness  value.  There  is  no  new  'genetic-material'  generated  in  the  selection 
process. 

Equation  (4.10)  shows  that  the  probability  of  a  design  being  selected  for 
participation  in  the  reproduction  process  depends  not  only  on  its  fitness,  but  also 
fitness  of  other  designs  in  the  population.  Furthermore,  the  value  returned  by  a 
fitness  function  is  not  always  an  "exact"  measure  of  fitness.  The  exact  value 
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returned  can  vary  greatly  depending  on  how  the  fitness  function  is  implemented. 
For  example,  let  a  population  in  a  minimization  problem  consist  of  four  designs 
which  have  objective  function  values  as  3,  5,  7,  and  9.  If  the  largest  value  9  is 
chosen  as  F\^^  (Eqn.  4.9),  then  the  fitness,  fj,  for  these  four  designs  will  be  6,  4,  2, 
and  0,  respectively.  The  probability  of  selection  of  these  four  designs  are  6/12, 
4/12,  2/12,  and  0/12,  respectively  (Eqn.  4.10).  The  design  with  the  lowest  fitness 
will  not  be  selected  for  participation  in  the  reproduction  process.  The  probability 
of  the  design  with  the  best  fitness  is  3  times  that  of  the  design  with  the  second 
worst  fitness.  Therefore,  a  higher  fitness  average  is  likely  to  occur  after  the 
reproduction  process.  This  gives  genetic  algorithms  power  to  guide  the  search 
towards  more  promising  design  areas.  However,  if  the  above  designs  have  the 
following  objective  functions  3,  5,  7,  and  85.  The  outstanding  high  objective 
function,  85,  could  result  from  multiple  causes,  for  instance,  it  may  be  due  to  a 
severely  violated  constraint.  If  the  value  85  were  selected  as  F*„,3^,  then,  the  fitness 
for  these  four  designs  will  be  82,  80,  78,  and  0.  In  each  selection,  the  probabilities 
of  the  four  designs  for  participation  in  the  reproduction  process  would  be  82/240, 
80/240,  78/240,  and  0/240.  Almost  no  differences  exist  between  the  designs  when 
considering  the  probabilities  of  participation  with  the  exception  of  the  worst  case. 
This  phenomenon  makes  it  difficult  for  good  designs  to  stand  out  from  the 
population,  and  consequently  the  genetic  search  would  not  concentrate  in  the 
more  promising  areas.  This  illustrates  a  drawback  of  induced  bad  fitness  as  caused 
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by  violated  constraints,  and  the  advantage  of  using  the  modified  penalty  P,  which 
was  formulated  in  Equation  (4.8). 

Fitness  proportionate  reproduction  (Eqn.  4.10)  has  been  used  almost 
exclusively  for  allocating  reproductive  trials  of  genetic  algorithms.  Reproduction 
strategies  that  differed  from  the  traditional,  proportionate  reproduction  plan  were 
implemented  [75,76]  in  an  effort  to  obtain  improved  performance  of  genetic 
searches.  Baker  [75]  reported  experiments  where  reproductive  trials  were 
allocated  according  to  the  rank  of  individual  strings  in  the  population  rather  than 
by  individual  fitness  values  as  related  to  the  population  average.  The  rank-based 
allocation  of  reproductive  trials  slowed  convergence.  The  method  gave  a  slow,  but 
more  accurate  optimization.  Whitley  [76]  presented  new  evidence  and  arguments 
which  suggested  that  allocating  reproductive  trials  according  to  rank  can  be  used 
to  speed  up  a  genetic  search  if  appropriately  implemented. 

4.2.2  Crossover 

The  process  of  reproduction  assures  that  more  copies  of  dominant  or  fit 
designs  will  be  present  in  a  population.  The  crossover  process  allows  for  an 
exchange  of  design  characteristics,  among  members  of  the  population  pool,  with 
the  intent  of  improving  the  fitness  of  the  next  generation.  This  is  similar  to  the 
transfer  of  genetic  material  in  biological  reproduction  processes  as  facilitated  by 
DNA  and  RNA  strings.  Crossover  is  executed  by  selecting  chromosomal  strings  of 
two  mating  parents,  randomly  choosing  sites  on  the  strings,  and  swapping  strings 
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of  O's  and  I's  between  these  chosen  sites.  There  are  a  number  of  crossover 
operations,  such  as  a  one-point  crossover,  two-point  crossover,  multi-point 
crossover,  segmented  crossover,  uniform  crossover,  and  shuffle  crossover.  The 
two-point  crossover  operator  was  used  in  this  research.  An  illustration  of  the  two- 
point  crossover  process  between  mating  parents  represented  by  20  digit  binary 
strings  is  as  follows: 

Parentl  =  11010100100110011010 
Parent2  =  01001110001011001000 

Childl  =  11001110001010011010 
Child2  =  01010100100111001000 

The  crossover  sites  on  the  parent  strings  are  indicated  by  an  understrike. 
A  probability  of  crossover  p,  is  defined  to  determine  if  crossover  should  be 
implemented.  In  the  case  of  the  more  fundamental  one-point  crossover  process, 
only  one  randomly  generated  site  is  created,  and  strings  after  the  selected  site  are 
then  swapped  between  two  partners.  More  information  about  the  implementation, 
nature,  and  evaluation  of  different  crossover  operators  is  presented  in  [77]. 
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4.2.3  Mutation 

Mutation  safeguards  the  genetic  search  process  from  a  premature  loss  of 
valuable  genetic  material  during  reproduction  and  crossover.  The  process  of 
mutation  is  simply  to  choose  few  members  from  the  population  pool  according  to 
the  probability  of  mutation  p„„  and  to  switch  a  0  to  1  or  vice  versa  at  a  randomly 
selected  mutation  site  on  the  selected  string.  Most  genetic  search  methods  give  a 
constant  probability  of  mutation,  p^.  A  method  employing  variable  probability  of 
mutation  in  the  genetic  search  was  proposed  by  Fogarty  [78].  The  study  measured 
the  performance  of  genetic  searches  in  which  the  p,„'s  were  a)  exponentially 
decreased  with  successive  generations,  b)  increased  over  the  bit  representation  of 
each  integer,  and  c)  was  a  combination  of  the  increasing  and  decreasing 
probability  methods.  The  variable  p„,  methods  were  compared  with  the 
performance  of  a  genetic  search  with  a  fixed  p„.  The  varied  probability  p^  was 
shown  to  significantly  improve  the  performance  of  the  genetic  search,  but  only  for 
cases  in  which  the  initial  populations  contained  all  O's  or  all  I's.  Randomly 
generated  initial  population  cases  showed  no  improvement  in  performance. 

Running  a  genetic  algorithm  entails  setting  a  number  of  initial  parameter 
values,  such  as  the  size  of  population,  probabilities  of  crossover,  p^,  mutation 
probability,  p^,  etc.  Determining  initial  conditions  which  suit  one's  problem  is  not 
a  trivial  task.  If  poor  initial  settings  are  used,  the  performance  of  a  genetic  search 
can  be  severely  impacted.  Two  useful  strategies  for  determining  effective 
crossover  and  mutation  probabilities  have  been  developed  [79,80].  A  strategy  to 
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adapt  the  probabilities  of  crossover  and  mutation  during  the  course  of  the  genetic 
search  was  implemented  by  Davis  [81].  This  method  showed  Hmited  success,  and 
more  efforts  were  needed  to  make  it  robust  and  effective.  The  performance  of 
genetic  search  methods  in  obtaining  the  optimum  ply  stacking  sequence  of  a 
laminated  composite  plate  under  a  concentrated  load  was  examined  for  a  number 
of  combinations  of  p^  and  p„,.  This  experiment  is  detailed  in  Chapter  5. 

4.2.4  Implicit  Parallelism 

It  is  important  to  emphasize  that  genetic  algorithms  use  probabilistic 
transition  rules  only  to  guide  a  highly  exploitative  search.  To  this  extent,  they 
should  not  be  considered  a  variant  of  the  random  walk  approach.  The  implicit 
parallelism  available  in  these  methods  is  significant  from  a  computational 
standpoint,  and  has  been  explored  [82].  If  a  binary  code  is  used  in  an  n-digit 
string  to  represent  a  design,  then  a  total  of  2"  variations  of  the  design  are 
available  in  the  representation  scheme.  However,  the  exploitative  power  of  OA's 
extends  beyond  this  if  one  takes  into  consideration  the  fact  that  the  presence  of  a 
0  or  1  at  some  key  location  along  the  string  may  be  significant  in  itself.  As  an 
example,  consider  the  design  of  a  two-bar  truss  structure  for  minimum  weight, 
with  each  cross-sectional  area  mapped  into  a  5-digit  binary  string.  The  maximum 
and  minimum  areas  would  be  denoted  as  follows. 
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=  11111 
=  00000 


If  the  bit  string  representation  of  the  design  is  obtained  by  stacking  the  binary 
representation  of  area  A,  after  Aj,  then  a  design  represented  by  1001000010 
indicates  a  value  of  Aj  close  to  its  maximum  and  A,  closer  to  its  minimum  value. 
Further,  if  we  use  a  '*'  to  indicate  that  a  location  along  the  string  can  occupy 
either  a  1  or  a  0,  then  a  string  of  the  form  1***0000**  conveys  the  same 
information.  In  this  manner,  one  can  think  of  the  following  schemata: 


J_J2     _  «Q«|***«|« 

H2  =  1***00***0 
H3  =  1*01***0** 


Our  design  of  areas  A^  and  Aj  could  belong  to  each  one  of  these  schemata.  In 
such  a  10-digit  representation,  there  are  a  total  of  3^°  (59049)  schemata  among  2'° 
(1024)  possible  unique  designs.  A  decrease  in  the  number  of  *'s  in  a  schema 
makes  it  more  specific.  In  this  context  one  defines  the  order  of  a  schema  denoted 
as  0(H)  as  the  number  of  O's  or  I's  in  a  schema.  The  order  of  Hi  is  3  and  that 
of  H2  and  H3  is  4.  Another  item  of  interest  is  the  defining  length  of  a  schema 
d(H),  which  is  taken  as  the  distance  between  the  first  and  last  specific  digit  on  the 
string;  d(Hl)  =  7,  while  d(H2)  =  9.  Schemata  with  higher  defining  lengths  are  more 
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likely  to  be  disrupted  during  crossover.  If  we  denote  the  number  of  schemata  of 
type  H  at  a  given  generation  t  as  m(H,t),  then,  due  to  reproduction  alone,  at 
generation  t+ 1  this  number  would  be  [82], 

f a\' 

where  f(H)  is  the  fitness  of  schema  H  and     is  the  average  fitness  of  the 
population.  Clearly,  schema  with  a  fitness  higher  than  the  average  fitness  of  the 
population  are  increased  exponentially  in  successive  generations.  The  affect  of 
crossover  and  mutation  is  disruptive.  At  the  same  time,  however,  these  processes 
contribute  to  the  evolution  of  new  characteristics  in  the  population.  In  terms  of 
the  probabilities  of  crossover  and  mutation,  p^  and  p^,  respectively,  and  the  string 
length  L,  the  growth  or  decay  of  a  schema  in  a  population  is  expressed  as  follows: 

m{H,t^l)  >  miH,t)l^(l-  p^^-  0{H)p\  (4.12) 


f 


L-1 


Here,  the  predicted  growth  due  to  reproduction  alone  has  been  modified  by 
considering  the  probability  of  surviving  the  crossover  and  mutation  transforms.  As 
IS  dear  from  this  expression,  low-order,  shorter  defining-length  schemata,  which 
are  also  fit,  increase  rapidly  in  the  population.  Holland  [34]  estimated  that  the 
actual  processing  of  n  strings  corresponding  to  the  population  size  results  in  the 
processing  of  O(n^)  schemata,  giving  the  approach  a  tremendous  computational 
advantage. 
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The  choice  of  population  size,  n,  for  a  genetic  search  represents  a 
fundamental  decision  in  GA  implementation.  A  population  with  too  few  members 
can  result  in  a  premature  convergence  of  the  genetic  search,  where  not  enough 
genetic  combinations  are  processed  to  converge  to  the  ideal  solution.  Premature 
convergence  resulted  from  an  effect  referred  to  as  genetic  drift  [79,83]  or  a 
preferential  selection  pressure.  An  overly  large  population  size  results  in  increased 
number  of  function  evaluations  before  significant  improvement  is  observed  in  the 
serial  evaluation  of  individuals.  Goldberg  [84,85]  provides  strategies  for 
determining  the  optimal  initial  population  sizes  of  binary-coded  serial  and  parallel 
genetic  algorithms.  The  performances  of  genetic  search  methods  were  computed 
for  a  number  of  structural  and  mechanical  optimal  design  problems  using 
different  population  sizes.  The  results  are  presented  in  Chapter  5. 

Premature  convergence  can  be  defined  as  an  event  where  a  population  of 
designs  achieves  a  high  level  of  uniformity  at  all  loci,  without  containing  near- 
optimal  design  characteristics.  This  undesirable  convergence  situation  usually 
results  from  stochastic  errors  associated  with  a  small  population  size.  Strings 
representing  important  schemata  are  not  usually  represented  sufficiently  when  a 
small  population  size  is  selected.  Finite  populations  eventually  converge  to  non- 
optimal  solutions  without  the  presence  of  selective  pressure.  This  phenomenon  is 
given  the  name  genetic  drift.  A  large  population  size  or  a  high  probability  of 
mutation  decreases  the  occurrence  of  genetic  drift.  Detailed  descriptions  of 
genetic  drift  and  of  the  Markov  analysis  of  genetic  drift  are  presented  in  [83]. 
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4.3  Representation  Schemes  for  Design  Space 
4.3.1  Continuous  Design  Variables 

An  m-digit  binary  number  representation  of  a  continuous  design  variable 
allows  for  2""  distinct  variations  of  that  design  variable  to  be  considered.  If  the 
description  of  a  design  variable  is  required  to  a  precision  of  A^,,  then  the  number 
of  digits  in  the  binary  string  may  be  estimated  from  the  following  relationship: 


2"'  > 


(4.13) 


Here,  Xl  and  x^^  are  the  lower  and  upper  bounds  of  a  continuous  variable  x.  As 
an  example,  if  one  considers  representing  a  variable  x  to  within  a  precision  0.1, 
and  the  lower  and  upper  bounds  on  x  are  0.0  and  0.7,  then  2'"  >  8.  This  gives  the 
result  that  m  =  3,  and  indeed  the  eight  3-digit  combinations  of  0  and  1,  enlisted  as 
000,  001,  010,      Ill,  can  be  assigned  to  numbers  0.0,  0.1,...0.7.  Note  that  an 
integer  value  of  m  will  not  yield  a  precision  of  0.14,  and  a  3-digit  string  would  still 
be  necessary  to  represent  the  six  numbers  0.0,  0.14,  0.28,     0.7.  There  would  be 
two  excessive  binary  representations  which  would  be  digested  by  distributing  them 
to  partial  discrete  design  space  as  described  in  later  sections.  It  is  important  to 
recognize  that  even  when  dealing  with  continuous  variables,  GA's  work  on  a 
discrete  representative  set  of  those  variables;  the  method  is  therefore  ideally 
suited  for  applications  to  problems  with  a  mix  of  continuous,  integer,  and  discrete 
variables. 
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The  binary  representation  of  population  characteristics  has  dominated  GA 
research.  Advantages  of  binary  representation  include:  analytical  simplicity  of 
binary  vectors,  the  elegance  of  genetic  operators  in  binary  strings,  requirements 
for  computational  speed,  and  most  importantly,  the  theoretically  motivated  desires 
to  maximize  the  number  of  schemata  sampled  by  a  given  set  of  individuals  in  a 
population  [34,82].  Binary  representation  is  not  an  exclusive  model  for  genetic 
algorithms,  and  higher-level  representations  in  both  theory  and  applications  are 
v^'orthy  of  investigation  [86]. 

4.3.2  Integer  Design  Variables 

Due  to  the  discrete  nature  of  the  binary  representation  schemes,  integer 
design  variables  can  be  simply  regarded  as  continuous  design  variables  with  a 
fixed  accuracy  A,  equal  to  1.  If  m  could  be  found  to  meet  (x,j  -  Xl)  =  2"'-l,  a 
one-to-one  correspondence  could  be  readily  established.  In  most  cases,  however, 
this  is  not  possible,  and  the  excessive  binary  strings  must  be  assigned  in  an 
appropriate  manner.  There  are  a  number  of  ways  in  which  this  is  done,  and  these 
are  described  as  follows. 

(i)  Penalty  approach:  In  this  approach,  the  smallest  number  m  which  meets  the 
inequahty 

2'"  >  {X,  -xj  .  1  (4.14) 
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is  computed.  Of  the  2"'  possible  m-digit  binary  strings,  a  unique  string  is  assigned 
to  each  of  the  N  integer  variables.  The  remaining  (2""  -  N)  strings  are  assigned  to 
out-of-bound  integers.  As  an  example,  in  representing  six  integers  between  0  and 
5,  the  computed  value  of  m  satisfying  the  above  inequality  is  3.  The  remaining 
two  binary  strings  are  assigned  to  out-of-bound  variables  as  follows, 

[0,  1,  2,  3,  4,  5,  6*,  7*] 

[000,  001,  010,  Oil,  100,  101,  110*,  111*] 

where,  '*'  indicates  an  out-of-bound  variable.  A  penalty  measure  is  then  allocated 
to  the  fitness  function  of  a  design  which  includes  the  out-of-bound  value  of  an 
integer  variable.  While  this  approach  yields  a  one-to-one  correspondence 
between  the  integer  variables  and  their  binary  representations,  careful 
consideration  must  be  given  to  the  magnitude  of  penalty  assigned  to  the  fitness 
function  due  to  the  presence  of  an  out-of-bound  variable.  Large  penalties  on  the 
fitness  will  adversely  affect  the  genetic  search,  making  it  difficult  to  distinguish 
between  good  and  average  designs. 

(ii).  Excessive-Distribution  Method:  In  this  approach,  m  is  first  computed  on  the 
basis  of  Eqn.  (4.14).  The  excessive  binary  representations  are  then  assigned  to 
integers  in  the  admissible  range,  whereby  one  or  more  integers  may  have  more 
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than  one  binary  representation.  For  the  case  described  above,  one  of  the  two 
binary  strings  could  be  assigned  to  integer  4  and  the  other  to  integer  5.  Although 
this  method  does  not  create  any  additional  design  constraints,  it  does  have  the 
effect  of  an  uneven  expansion  in  the  design  space.  The  affect  of  this  partial 
expansion  is  hard  to  predict  without  apriori  information  of  the  objective  function 
space.  Clearly,  the  expansion  results  in  a  larger  number  of  schemata  in  the  design 
space.  If  the  schemata  happen  to  belong  to  the  poorer  regions  of  the  design  space, 
the  convergence  of  genetic  search  would  be  adversely  influenced.  One  method  of 
avoiding  this  problem  is  to  distribute  the  excess  binary  strings  evenly  along  the 
extent  of  the  feasible  integer  space.  As  an  example,  if  twenty  excessive  binary 
strings  are  to  be  distributed  evenly  among  one  hundred  integers,  one  excess  binary 
representation  can  be  assigned  to  every  fifth  consecutive  integer. 

The  uneven  distribution  described  above  may  be  avoided  by  creating  excess 
strings  on  purpose  to  obtain  a  similar  number  of  binary  representations  for  each 
integer.  For  the  example  described  above,  two  integers  (33%)  have  twice  the 
number  of  binary  representations  (100%  more)  than  the  other  four  (67%)  of  the 
integers.  If  m  =  4,  we  would  have  a  total  of  16  binary  strings  to  be  assigned  to  six 
integers.  In  this  case,  each  integer  could  be  assigned  2  strings,  and  four  integers 
(67%)  could  be  assigned  one  additional  string  (50%  more)  than  the  remaining 
two  (33%)  integers.  With  increasing  value  of  m,  the  disparity  in  distribution  can 
be  removed  at  the  price  of  increasing  the  number  of  schemata  to  be  explored. 
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4.3.3  Discrete  Design  Variables 

Discrete  type  design  variables  are  characterized  by  an  uneven  spacing 
between  two  consecutive  values.  With  the  unique  nature  of  genetic  search, 
wherein  gradient  information  is  not  required,  these  variables  can  be  handled  in 
the  same  manner  as  continuous  or  integer  variables.  Mapping  of  these  variables 
is  a  two-stage  process.  In  the  first  stage,  discrete  variables  are  mapped  to  an 
equivalent  number  of  integer  variables.  Then,  techniques  of  mapping  integer 
variables  into  binary  strings  described  in  the  preceding  sections  are  applicable 
with  no  additional  manipulation.  An  illustration  of  this  process  is  the  mapping  of 
a  set  of  eight  discrete  variables  into  binary  strings. 

[2.4,  3.76,  5.96,  8.25,  9.37,  13.70,  20.55,  24.0] 
[  1     2      3      4      5       6       7       8  ] 
[000  001    010   Oil    100     101     110    111  ] 

A  design  space  with  a  mix  of  continuous,  discrete,  and  integer  variables  can  be 
represented  as  required  in  genetic  search  by  connecting,  head-to-tail,  binary  string 
equivalents  of  these  variables  as  described  above. 

In  genetic  search,  the  average  fitness  of  a  population  is  increased  over 
generations  of  evolution.  In  an  unconstrained  maximization,  the  objective 
function  can  be  chosen  as  the  fitness  function.  If  an  unconstrained  minimum  is  the 
objective,  the  fitness  must  be  revised  to  be  the  inverse  of  the  objective  function. 
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As  stated  earlier,  another  alternative  is  to  choose  the  maximum  objective  function 
value  of  all  members  of  the  population,  and  subtract  individual  objective  function 
values  from  this  quantity  to  obtain  the  corresponding  fitness  values.  In 
constrained  optimization,  designs  with  violated  constraints  are  considered  less  fit. 
Hence  in  constrained  minimization,  a  penalty  corresponding  to  the  constraint 
violations  is  appended  to  the  objective  function  much  as  in  the  exterior  penalty 
function  approach.  In  the  present  work,  the  traditional  penalty  function  approach 
is  compared  to  a  revised  formulation  discussed  earlier,  which  is  characterized  by 
the  imposition  of  an  upper  bound  on  the  penalty.  Without  such  a  bound,  a  highly 
penalized  infeasible  design  and  a  good  feasible  design  would  be  difficult  to 
separate,  and  the  genetic  search  would  deteriorate  into  a  random  search. 

4-3.4  Performance  Evaluation  of  Representation  Schemes 

DeJong  [79]  defines  online  and  offline  performance  measures  to  gage  the 
efficiency  of  a  genetic  search.  While  online  performance  is  an  average  of  all 
evaluations  of  fitness,  and  is  therefore  indicative  of  how  well  the  entire  population 
adapts,  the  offline  parameter  is  an  average  of  the  current  best  evaluations  of 
fitness.  These  two  parameters,  in  addition  to  the  average  fitness  of  a  population 
at  each  generation  of  evolution,  were  used  to  measure  the  performance  of  genetic 
search  with  the  different  representation  schemes  described  in  preceding  sections. 

A  test  problem  for  this  purpose  was  chosen  as  follows: 

Minimize      F^(x)  =  x'^  +  jc,'  +  x,"  (4.15) 
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where  Xj,  x,,  and  X3  could  assume  continuous  values  ranging  between  -512  and  511. 
A  one-to-one  correspondence  can  be  obtained  by  using  a  10-digit  binary  string  (2'° 
distinct  representations).  String  lengths  of  12,  14,  16,  18  and  20  were  also  used, 
yielding  4,  16,  64,  256  and  1024  binary  representations  per  integer  variable, 
respectively.  A  fixed  population  size  of  50  was  used  in  conjunction  with 
probabilities  of  crossover  and  mutation  set  to  0.6  and  0.005,  respectively.  Both 
online  and  offline  performance  of  each  of  these  cases  was  obtained.  The  results 
show  improvement  in  both  performance  measures  with  an  increase  in  the  number 
of  representations  for  each  admissible  value  of  the  integer  variables  (20-digit 
string  performed  the  best).  This  was  attributed  to  two  possible  factors. 

a.  More  schemata  are  available  in  such  larger  string  representations,  thereby 
increasing  the  probability  of  having  a  better  initial  population. 

b.  In  populations  where  there  is  exactly  one  representation  for  each 
admissible  value  of  the  integer  variable,  there  is  the  risk  of  losing  that 
variable  very  easily  during  genetic  transformations.  Note  that  each  0  or  1 
along  the  string  very  specifically  denotes  an  integer  variable,  and  is 
susceptible  to  elimination  by  even  a  simple  mutation  transform.  These 
advantages  must  be  weighed  against  the  slower  improvement  in  the 
population,  a  direct  consequence  of  schema  disruption  with  longer  binary 
strings. 
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A  modification  of  this  problem  was  a  change  in  the  design  variable  upper 
bound  to  238.  For  string  length  of  10,  there  were  excess  binary  representations, 
which  were  treated  by  the  methods  discussed  earlier.  In  the  penalty  approach,  a 
penalty  (O(IO^))  was  imposed  on  any  design  containing  design  variable  values 
between  239  and  511.  This  choice  of  penalty  parameter  is  not  advised  because 
high  penalties  on  the  objective  function  resulted  from  the  design  variable  range 
violations  will  slow  the  convergence  of  the  genetic  search.  The  magnitude  of  the 
penalty  for  the  design  variable  range  violation  should  be  appropriately  selected  in 
accordance  to  the  average  magnitude  of  fitness  functions.  An  arbitrary  assignment 
of  penalty  parameters  should  be  avoided.  Another  genetic  search  with  an  initial 
penalty  (0(10"*)),  and  which  was  linearly  decreased  in  every  generation,  was  also 
used  with  more  satisfactory  results. 

In  the  excess  distribution  method,  two  schemes  were  implemented.  In  the 
first,  the  273  excess  representations  were  distributed  evenly  over  the  range  -512  « 
238;  in  the  second,  they  were  all  distributed  to  integers  on  the  side  of  238.  Four 
other  experiments  were  conducted,  with  string  lengths  of  11,  12,  13  and  14,  and  all 
excessive  representations  were  distributed  evenly  for  the  entire  range  of  integer 
variables.  It  can  be  shown  that  a  string  length  of  13  yields  the  most  even 
distribution  with  91%  of  the  integer  variables  assigned  only  about  10%  more 
representations  than  the  remaining  9%  of  the  integer  variables. 

Some  general  conclusions  emerged  from  these  numerical  experiments. 
Although  conventional  wisdom  in  genetic  search  advocates  the  use  of  smaller 
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length  strings,  for  a  discrete  design  space  it  appears  advantageous  to  have  more 
than  a  one-to-one  correspondence  between  the  integer/discrete  design  variables 
and  their  binary  representations.  If  a  one-to-one  correspondence  is  chosen  with 
only  a  few  excess  representations,  then  the  design  variable  representation  can  be 
based  on  either  the  penalty  concept  or  the  excessive  distribution  approach.  The 
penalty  approach  is  highly  sensitive  to  the  choice  of  the  penalty  parameter,  and  is 
generally  not  recommended  unless  information  about  the  magnitude  of  the 
objective  function  is  available.  Of  the  excessive  distribution  method,  the  even 
distribution  schemes  perform  better  than  the  one-sided  distribution  methods.  This 
conclusion  would  not  be  valid  only  if  excessive  distribution  is  done  in  those 
regions  of  the  design  space  where  the  optimal  point  is  located.  When  an 
excessive  distribution  method  is  used  for  a  design  variable,  the  design  space 
corresponding  to  this  variable  will  be  partially  expanded.  With  larger  length 
binary  numbers,  the  excessive  strings  can  be  more  evenly  distributed  and  result  in 
an  even  expansion  over  all  variables.  This  approach  also  makes  available  a  larger 
number  of  schemata  for  each  design  variable. 


4.4  Multicriterion  Design 
A  relatively  simple  approach  to  account  for  the  multicriterion  nature  of  the 
problem  in  optimal  design  is  through  the  utility  function  method  [74].  A  utility 
function  is  defined  by  assigning  a  relative  importance  to  an  individual  criterion  in 
relation  to  other  criteria.  For  an  m-criterion  optimization  problem,  a  practical 
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definition  of  the  utility  function  can  be  shown  to  be, 


U  =  YlW^.(x)  (4.16) 


1=1 


where      serves  as  the  weighting  factor  for  each  objective  function  Fj(x).  Since 
the  magnitude  of  each  objective  criterion  can  be  quite  different,  a  scaling  strategy 
must  be  implemented.  Here,  the  following  scaled  form  of  the  utility  function  was 
used, 


-         m  p 

U  =  Y:  W~  (4-17) 

'=1  f: 


where  F  '  are  the  scaling  parameters  for  each  candidate  criterion. 

With  this  definition,  the  multicriterion  optimal  design  problem  is  changed 
into  a  single  objective  optimal  design  problem,  and  can  be  solved  in  a  relatively 
straightforward  manner.  In  a  situation  where  the  multicriterion  optimal  design  is 
required  for  a  number  of  different  weighting  combinations,  a  genetic  search  would 
have  to  be  performed  for  each  of  these  cases  separately.  One  of  the  principal 
contributions  of  the  present  research  is  an  exploration  of  methods  wherein  single 
runs  of  genetic  search  can  be  used  to  simultaneously  generate  a  family  of  Pareto 
optimal  designs  corresponding  to  different  weighting  coefficients  of  individual 
candidate  criterion.  Two  distinct  strategies  of  achieving  this  objective  are 
proposed  next. 
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4.4.1  Sharing  Function  Approach 

Goldberg  and  Richardson  [87]  introduced  the  concept  of  sharing  functions 
into  genetic  algorithms  to  simultaneously  locate  relative  optima  in  a  multimodal 
design  space.  As  the  name  suggests,  the  approach  is  based  on  a  concept  of  shared 
resources  among  distinct  sets  of  population.  Each  such  set  converges  to  one 
relative  optimum,  and  in  doing  so,  maximizes  its  payoff.  The  principle  of  sharing 
is  implemented  by  degrading  the  fitness  of  each  design  in  proportion  to  the 
number  of  designs  located  in  its  neighborhood  through  the  use  of  sharing 
functions.  The  extent  of  sharing  is  controlled  by  a  sharing  parameter  o^^^,  in  terms 
of  which  the  sharing  function  is  defined  as  follows: 


sh  ^ 


0 


4ij  <  o 


sh 


(4.18) 


otherwise 


Here  djj  is  a  metric  indicative  of  the  distance  between  designs  i  and  j.  If  two 
designs  i  and  j  are  very  close,  djj  is  almost  zero,  and  0(d,j)  «  1.  If  d^  >  a.^,  <p{d,) 
=  0.  The  fitness  of  a  design  i  is  modified  as, 


I  m 


(4.19) 


7=1 
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where  m  is  the  number  of  designs  located  in  vicinity  of  the  i""  design.  Note  that 
the  raw  fitness  value  of  each  design  is  decreased  due  to  the  presence  of  a  large 
number  of  designs  in  its  vicinity. 

If  the  distance  metric  d^  is  evaluated  in  the  decoded  design  space,  the 
sharing  is  called  a  phenotypic  sharing.  The  distance  metric  djj  can  be  computed  as 
follows: 


In  order  to  adapt  the  sharing  principle  into  an  approach  where  genetic 
search  can  simultaneously  locate  the  optima  corresponding  to  different  weighting 
combinations,  the  weighting  variables  were  included  in  the  set  of  design  variables. 
As  in  mathematical  programming  based  multicriterion  design,  the  criterion 
weights  W.  are  required  to  meet  the  following  requirement. 

m 

E      =  1  (4.21) 

(=1 


As  an  illustration,  if  9  different  weighting  combinations  [Wj  =  0.1,  W,  =  0.9], 
[Wi  =  0.2,  W2  =  0.8],      [Wi  =  0.9,  W,  =  0.1]  are  considered  for  a  two-criterion 
optimal  design  problem,  the  weighting  variable  can  be  defined  as  follows: 
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xw     =  (1,2,  ,  9) 

=  ([0.1,0.9],  [0.2,0.8],  ...[0.9,0.1]) 

This  weighting  variable      is  then  regarded  as  one  additional  design 
variable  taking  on  values  between  1  and  9,  with  a  step  size  of  1.  This  variable  is 
of  integer  type,  and  representation  schemes  discussed  earlier  can  be  used  in  a 
routine  manner. 

For  the  above  two-criterion  problem,  one  weighting  design  variable,  Xy^,, 
was  needed  to  represent  the  nine  weighting  sets.  If  the  number  of  criteria  is 
greater  than  two,  then  two  approaches  can  be  taken  to  represent  the  weighting 
variable  sets.  The  first  approach  introduces  m-1  integer  variables,  x^,  into  the 
design  variable  set  of  an  m-criterion  problem.  Consider  the  following  nine 
different  weighting  combinations  for  a  three-criterion  problem: 

[Wi=0.2,  \V2  =  0.3,  W3  =  0.5],  [Wj  =  0.2,  W,  =  0.4,  W3  =  0.4], 

[W,=0.2,  W.-0.5,  W3  =  0.3],  [Wi  =  0.3,  W2  =  0.3,  W3  =  0.4], 

[Wj  =  0.3,  W,  =  0.4,  W3  =  0.3],  [W,  =  0.3,  W,  =  0.5,  W3  =  0.2], 

[Wj  =  0.4,  W,  =  0.3,  W3  =  0.3],  [Wi  =  0.4,  W,  =  0.4,  W3  =  0.2], 
[W,=0.4,  W2  =  0.5,  W3  =  0.1]. 

Two  weighting  variables  x^i  and  x^,  were  used  to  represent  W,  and  W„ 
respectively.  The  value  of  W3  in  each  weighting  set  was  then  determined  by  Eqn. 


90 

(4.21).  The  second  approach  utilizes  a  single  weighting  variable  which  represents 
the  possible  combinations  of  the  optimization  problem.  Weighting  combinations 
are  arranged  in  a  series,  and  the  weighting  variable  components  map  to  each 
combination  of  the  problem  criteria.  For  the  above  nine-criterion  example, 
the  following  variable  map  was  generated: 

xw  =  (1,  2,  9) 

=  ([0.2,0.3,0.5],  [0.2,0.4,0.4],  [0.4,0.5,0.1]) 

The  binary  coded  weighting  variable  and  binary  representations  of  true 
design  variables  are  connected  head-to-tail  to  form  a  binary  string  representing  a 
design  with  a  specific  weighting  combination.  In  order  to  generate  optimal  designs 
corresponding  to  all  weighting  combinations,  the  sharing  principle  is  applied 
selectively  on  the  weighting  variable.  The  distance  metric  d,j  between  two  designs 
Xj  =  [x\v,i,  Xj  j,  x,  I,     Xp  j]  and  Xj  =  [x^j,  x^j,  X2j,     XpJ  is  based  only  on  the 
variable  representing  weight  combinations  as  follows. 

With  appropriate  choices  of  sharing  parameter  o^^  and  population  size,  sharing 
directs  the  genetic  search  toward  locating  optimal  designs  corresponding  to  most 
of  weighting  combinations.  In  order  both  to  speed  convergence  and  to  impart 
stability  to  the  genetic  search,  it  was  found  to  be  necessary  to  introduce 
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restrictions  in  mating  between  members  belonging  to  different  sub-sets  of  the 
population  [88].  A  parameter  of  mating  restriction,  a^^,,  was  defined.  This 
confined  the  mating  process  among  those  members  for  which  the  metric  distance 
as  measured  along  the  weighting  variables  was  within  the  radius  a^^,.  Some 
caution  is  necessary  in  this  approach  as  it  is  highly  reminiscent  of  in-breeding,  and 
sources  of  diversity  in  the  population  must  be  introduced.  It  is  important  to 
include  in  future  work,  mating  restrictions  other  than  those  based  entirely  on  the 
weighting  coefficient  combinations.  Since  the  population  of  designs  in  each 
specific  weighting  group  is  determined  on  the  basis  of  the  penalized  fitness  of  the 
entire  population,  a  better  average  fitness  for  a  specific  weighting  group  will 
eventually  result  in  that  group  attracting  more  members  in  successive  generations. 
The  fitness  or  the  overall  utility  function  has  to  be  carefully  scaled  so  as  to 
distribute  designs  evenly  over  each  weighting  combination,  and  for  each  weighting 
group  to  possess  roughly  similar  exploitative  abilities  in  the  genetic  search. 

4.4.2  Vector-Evaluated  Approach 

This  approach  has  its  parallel  in  a  method  proposed  by  Schaffer  [89], 
referred  to  as  the  vector  evaluated  genetic  algorithm  (VEGA)  approach,  and 
intended  for  multiobjective  optimization.  Given  the  fuced  population  size  as  P 
and  the  number  of  optimization  criteria  as  N,  each  member  of  the  population  was 
evaluated  for  its  fitness  as  it  relates  to  each  of  the  N  criteria.  For  each  of  the  N 
criteria,  subpopulations  of  size  P/N  were  extracted  from  the  pool  on  the  basis  of 
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the  fitness  function  for  the  criterion  under  consideration.  This  sorting  was  done 
in  a  rotational  sequence,  where  one  member  was  selected  for  each  criterion, 
followed  by  a  second  round  of  selection.  The  selection  was  performed  using  the 
traditional  genetic  selection  (reproduction)  method,  which  ensured  the 
preservation  of  individual  designs  which  performed  above  average  on  any 
criterion.  Subpopulations  were  then  shuffled  and  grouped  into  a  population  which 
was  subjected  to  the  usual  genetic  transfer  processes. 

The  method  proposed  in  the  present  research  for  simultaneously  generating 
a  family  of  Pareto  optimal  designs  draws  upon  such  a  vector  evaluated  approach, 
and  is  summarized  as  follows.  A  utility  function  is  first  defined  to  obtain  a  scalar 
measure  of  fitness  for  a  given  set  of  weighting  parameters  W,""  as 

m 

U,  =  E       F,  (4.23) 

where  the  superscript  k  denotes  the  k-th  weight  combination.  Each  set  of  weight 
combination  for  which  an  optimal  design  is  needed  constitutes  a  component  in  the 
N-vector  of  utility  functions.  The  entire  population  of  P  designs  was  evaluated  for 
each  of  the  weight  combinations.  A  subpopulation  of  size  P/N  for  each  weighting 
combination  was  formed  as  described  above.  Two  distinct  strategies  were 
adopted  in  this  context.  In  the  first  approach,  designs  were  selected  from  the  pool 
as  those  that  were  deemed  most  fit  for  that  weight  combination.  Selection  of  a 
design  for  one  subpopulation  did  not  exclude  its  subsequent  selection  in  another 
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subpopulation.  In  the  second  approach,  an  exclusionary  principle  was  applied 
wherein  a  design  could  contribute  its  genetic  makeup  to  only  one  subpopulation 
and  was  not  considered  for  mating  in  another  subpopulation.  Both  of  these 
processes  result  in  N  subpopulations,  each  containing  P/N  designs  with  superior 
measures  of  fitness  corresponding  to  a  specific  weight  combination.  For  both 
methods,  in  general,  this  step  simply  represents  a  regrouping  with  no  genetic 
transformations  invoked  on  the  designs.  In  particular,  in  the  exclusionary  scheme, 
the  entire  genetic  pool  is  preserved  and  simply  divided  up  in  the  most 
opportunistic  manner  to  contribute  towards  evolution  of  designs  for  specific 
requirements.  After  the  selection  of  each  subpopulation  is  completed,  genetic 
transformations  such  as  reproduction,  crossover  and  mutation  are  invoked 
separately  in  each  subpopulation.  The  process  of  population  classification  is 
repeated  for  the  next  generation  in  the  same  manner.  In  the  discussion  of  sample 
results,  these  two  strategies  are  denoted  as  VEj^,  and  VE^^„  respectively. 

Due  to  the  possible  small  size  of  each  subpopulation,  genetic  drift  might 
result,  with  associated  undesirable  characteristics  such  as  convergence  to  a 
sub-optimal  design.  A  possible  approach  to  prevent  such  a  premature 
convergence  is  to  introduce  a  more  moderate  sized  population  for  genetic 
evolution.  The  underlying  assumption  here  is  that  a  larger  population  is  likely  to 
be  more  diverse,  and  hence  offers  a  higher  possibility  of  locating  the  true 
optimum.  This  can  be  achieved  in  a  strategy  where  the  crossover  transformation 
is  a  two-stage  process.  In  the  first  stage,  a  crossover  is  performed  between  a 


94 

subset  of  members  within  the  same  subpopulation  (akin  to  inbreeding).  The 
remaining  members  from  each  subset  forms  a  larger  gene  pool  where  a  second 
stage  of  crossover  is  conducted.  In  this  second  stage,  a  mating  restriction  must  be 
applied  on  the  variable  denoting  the  weighting  combination  associated  with  a 
design,  so  as  to  prevent  crossovers  between  designs  corresponding  to  radically 
different  weight  combinations. 


CHAPTER  5 

NUMERICAL  IMPLEMENTATION  OF  GENETIC  ALGORITHMS 


This  chapter  describes  the  numerical  implementation  of  genetic  algorithms 
in  single  and  multicriterion  optimization  problems.  Several  single-criterion 
optimization  problems  in  structural  and  mechanical  system  design  are  presented 
in  the  first  part  of  the  chapter.  The  second  section  describes  a  number  of 
multicriterion  optimization  problems  in  representative  structural  systems.  These 
single  and  multicriterion  problems  were  formulated  such  that  they  contained  a 
mix  of  integer,  discrete  and  continuous  design  variables.  The  two  variations  of 
genetic  search  described  in  the  previous  chapter  were  used  to  simultaneously 
obtain  families  of  Pareto  optimal  solutions. 

5.1  Single-Criterion  Optimization 

5.1.1  Riveted  Lap  Joint 

The  first  problem  involved  the  design  of  a  lap  joint  between  two  steel 
plates  in  which  the  rivet  size  and  the  number  and  arrangement  of  the  rivet  pattern 
were  considered  as  design  variables.    The  configuration  of  the  plates  and  the 
rivets  is  shown  in  Figure  5.1.  The  number  of  rows  parallel  to  side  AB  is 
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Figure  5.1  Geometry  of  riveted  lap  joint. 
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represented  by  an  integer  variable  Xj  with  permissible  values  between  1  and  32. 
The  number  of  the  rivets  in  each  row  x,  was  an  integer  variable  and  was  allowed 
to  assume  values  between  1  and  128.  The  diameter  X3  of  all  rivets  was  assumed 
to  be  the  same,  and  is  chosen  from  a  conmiercially  available  set: 


X3  =  [6,8,10,12,14,16,18,20,22,24,27,30,33,36,40,45]  (mm) 


This  variable  is  therefore  of  the  discrete  type.  The  objective  of  this  optimization 
was  to  maximize  the  efficiency  of  the  joint,  defined  as  the  ratio  of  the  strength  of 
joint  to  the  strength  of  the  plate.  The  strength  of  the  joint  is  obtained  as  the 
minimum  of  the  shearing  failure  strength  P,,  tension  failure  strength  P„  and 
compression  or  bearing  failure  P^.  These  can  be  formulated  as  follows  [90,91]: 


2  T 

TTX-<X,X^  — 
J    1  2  ^ 


nx  2X^X2 


4[1.06+0.126(Xj-3)] 


Jfx,  <  3 


,otherwise 


(5.1) 


P,  =  (2000  -  x^^)ta^ 


(5.2) 


P.  = 


[1.06+0.26(Xj-3)] 


,z/xj  <  3 
^otherwise 


(5.3) 
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where. 


T  =  80  MPa  a,  =  90  MPa  =  120  MPa 

Stress  concentrations  due  to  the  close  placement  of  any  two  rivets  were  avoided 
by  the  imposition  of  the  following  linear  constraints  [90]: 

SxjXj  +       <  500 

(5.4) 

SxjX,  +       <  2000 


Further,  the  strength  of  the  plate  was  specified  as  2700  kN  in  the  present  analysis. 
For  the  genetic  search,  population  sizes  of  30  and  60  were  considered,  with  two 
distinct  penalty  function  formulations.  The  total  number  of  allowed  values  for 
design  variables  Xj,  x,,  and  X3  were  32  (2^),  128  (2''),  and  16  (2^),  respectively. 
Design  variables  Xj,  x,,  and  X3  were  then  represented  as  binary  strings  of  5,  7  and 
4,  digits,  respectively.  Results  of  the  genetic  search  simulation  using  a  maximum 
bound  on  the  constraint  penalty  are  shown  in  Figures  5.2  and  5.3;  they  are 
marked  with  an  asterisk  to  distinguish  them  from  those  obtained  with  an 
unbounded  penalty  formulation. 

A  population  size  of  30,  when  used  with  a  bounded  constraint  penalty, 
performed  most  favorably.  An  optimum  design  of  Xj*  =  5,  x,'  =  13,  and  X3'  =  27mm 
was  obtained  at  the  58th  generation,  and  resulted  in  a  joint  efficiency  of  82.45%. 
The  average  efficiency  of  the  joint  for  both  population  sizes,  and  using  no  bounds 
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Figure  5.2  Iteration  history  of  lap  joint  efficiency-online  and  offline  measures. 
(*  indicates  use  of  bounded  penalty). 
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on  the  constraint  penalty,  actually  decreased  over  generations  of  evolution.  This 
underscores  the  need  for  a  careful  selection  of  penalty  functions  in  the  presence 
of  constraints. 

5.1.2  Series  Stacked  Belleville  Spring 

A  second  example  deals  with  the  sizing  and  selection  of  the  number  of  the 
belleville  springs  to  be  stacked  in  series  to  maximize  the  load  carrying  capacity 
with  a  minimum  attendant  deformation.  A  Belleville  spring  is  shown  in  Figure 

5.4,  and  the  configuration  for  a  series-stacked  Belleville  spring  is  shown  in  Figure 

5.5.  The  compressive  stress  at  the  convex  inner-diameter  corner  would  be  the 
critical  consideration  in  such  a  statically  loaded  situation.  The  load  deflection 
relation  is  obtained  as  follows  [92,93], 


where  d  is  the  deflection  under  load  P.  The  expression  for  maximum  stress  is 
obtained  as, 


P  = 


-[ih-f){h-0.5f)t  .  P] 


(5.5) 


(l-v)Mfl- 


-[C,(h  -  0.5/)  .  C,r] 


(5.6) 


(1  -  v-)Ma- 


where. 


Figure  5.5  Series  stacked  Belleville  spring. 
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C,  = 


M 


6 

R 

-  1 

-nlnR 
6 

\R  -  1 

nlnR 

2 

6 

{R  -  1)2 

nlnR 

(5.7) 


(5.8) 


(5.9) 


a  =  0.5  O.D. 
R  =  O.D./I.D. 

f  =  vertical  deflection  of  the  spring 

H  =  h  +  t  cos  q,  the  height  of  a  single  spring 

N  =  number  of  springs  in  stack 

Hj  =  N*H,  total  height  of  the  stacked  springs 


In  the  above  equations,  a  value  of  Poisson's  ratio  of  v  =0.3  and  Young's  modulus 
E  =  2.06  MN/mm-  were  used.  It  is  important  to  note  that  the  load-deflection 
relations  are  only  valid  for  d  <  0.85h.  Consequently,  the  objective  of  the 
optimization  is  to  find  the  maximum  load  at  d  =  0.8h,  and  subject  to  other 
constraints.  When  an  h/t  ratio  greater  than  1.3  is  used  in  a  stack,  the  load 
deflection  relation  tends  to  be  erratic,  as  some  of  the  springs  have  a  tendency  to 
snap  through  the  flat  position.  To  avoid  this  problem,  a  constraint  of  h/t  <  1.3 
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was  imposed.  The  maximum  allowable  working  stress  was  set  to  1400  MPA  and  a 
restriction  on  R  to  be  at  least  1.5  was  also  included. 

The  design  variables  Xj,  x,,     X5  are  the  thickness  t,  outer  diameter  O.D., 
inner  diameter  I.D.,  number  of  springs  N,  and  disk  cone  angle  6,  respectively. 
The  first  three  were  treated  as  discrete  variables  selected  from  the  following  sets. 


0.  D.  s  [200,  201,  202,...263]  (mm) 

1.  D.  s  [112,  113,  114,...175]  (mm) 

t     =  [4.8,  5.5,  6.05,  6.5,  7.1,  7.49,  9.4,  11.25]  (mm) 


The  number  of  springs  was  an  integer  variable  ranging  from  ]  to  16;  the  cone 
angle  0  was  a  continuous  variable  bounded  between  1  deg  and  20  deg.  It  was 
represented  to  an  accuracy  of  0.1  in  the  binary  scheme.  Population  sizes  of  50  and 
100  were  considered  in  the  numerical  simulation,  with  constraints  treated  both  by 
a  general  penalty  and  a  bounded  penalty  formulation.  Binary  strings  of  lengths  3, 
6,  6,  4,  and  10  were  chosen  for  design  variables  Xj  through  x^.  The  choices  of 
binaiy  strings  of  lengths  3,  6,  6,  and  4  for  x^,  x^,  X3,  and  x,  resulted  in  one-to-one 
mapping  between  the  actual  design  variable  space  and  its  corresponding  binary 
representation  space.  The  choice  of  a  binary  string  of  length  10  resulted  in  an 
even  distribution  where  1024  binary  representations  were  assigned  to  200 
candidate  Xj  values.  This  meant  that  5  representations  were  almost  uniformly 
distributed  to  each  design  value.  The  results  obtained  from  the  simulation  are 
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shown  in  Figures  5.6  and  5.7,  where,  as  before,  results  corresponding  to  the 
bounded  penalty  formulation  are  denoted  with  an  asterisk.  This  bounded  penalty 
formulation  yielded  the  most  stable  genetic  search  when  working  with  a 
population  size  of  50.  At  the  317th  generation,  an  optimal  design 

X*  =  [9.4  (mm),  232  (mm),  130  (mm),  7,  8.1  (deg)] 

was  obtained,  corresponding  to  a  maximum  load  of  142,211  lbs. 

5.1.3  Eleven  Bar  Planar  Truss 

A  minimum  weight  truss  design  with  constraints  on  displacements 
comprised  the  third  test  problem.  The  structure  and  the  loading  are  shown  in 
Figure  5.8.  The  cross-sectional  areas  of  each  of  the  eleven  members  of  the  truss 
are  discrete  variables  chosen  from  the  following  available  sections: 

A'  =    [0.111,  0.141,  0.196,  0.250,  0.307,  0.391,  0.442,  0.563, 

0.602,  0.766,  0.785,  0.994,  1.000,  1.228,  1.266,  1.563]  (in-). 

Material  properties  of  steel  with  specific  weight  p=  0.283  Ib/in^  and  Young's 
modulus  E=  30x10^  psi  were  prescribed  for  this  application.  Binary  string  length 
of  four  was  used  to  represent  the  design  variables;  this  results  in  exactly  one 
binary  string  for  each  variable.  The  displacement  constraints  were  defined  as 
follows: 
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Figure  5.6  Iteration  history  of  maximum  spring  Ioading--online  and  offline  measures 
(*  indicates  use  of  bounded  penalty). 
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Figure  5.8  Eleven  bar  planar  truss. 
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Ui  <  0.5  in. 

u,  <  0.5  in.  (5.10) 
U3  <  0.2  in. 

The  traditional  penalty  function  formulation  was  used  with  population  sizes 
of  60  and  120;  for  the  bounded  penalty  function  approach,  sizes  of  50  and  100 
were  considered.  These  results  are  summarized  in  Figures  5.9  and  5.10, 
respectively.  As  in  previous  problems,  the  bounded  penalty  formulation  yielded 
better  performance.  An  optimal  weight  of  82.61  lbs  was  obtained  with  the 
corresponding  design  variable  vector  as  follows: 

X*  =    [0.442,  0.602,  0.307,  0.442,  0.602,  0.111,  0.111,  0.602,  0.602,  0.442]  (in") 

5.1.4  Laminated  Composite  Plate  Design 

A  twenty-ply  fiber-reinforced  laminated  composite  plate  subjected  to  a 
central  concentrated  load  was  chosen  as  the  fourth  test  problem.  The  optimum 
stacking  sequence  of  the  twenty  unidirectional  laminae  was  determined  to 
minimize  the  central  displacement.  A  simply  supported  laminated  composite 
plate  under  a  concentrated  force  is  shown  in  Figure  5.11.  Two  laminated 
composite  plates,  one  square  and  the  other  rectangular,  were  considered  in  this 
problem.  Each  side  of  the  square  plate  was  0.5m.  The  rectangular  plate  was 
0.6m  by  0.3m.  The  0.6m  side  of  the  rectangular  plate  was  in  the  x-direction. 
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Figure  5.9 


Iteration  history  of  structural  weight-online  and  offline  measures 
(*  indicates  use  of  bounded  penalty). 
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Figure  5.10  Iteration  history  of  structural  weight-average  value  (*  indicates  use 
bounded  penalty). 
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Figure  5.11 


Simply  supported  laminated  plate  under  concentrated  load. 
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which  defined  the  zero  degree  direction  of  the  fiber  orientation.  A  vertical 
concentrated  force  of  lOON  was  apphed  at  the  center  of  both  plates.  The 
laminated  plate  was  assumed  to  be  symmetric  with  respect  to  its  middle  plane. 
This  impHed  that  only  plies  above  the  plane  of  symmetry  needed  to  be  considered 
in  the  design  variable  set.  These  ten  plies  were  further  assumed  to  consist  of  two 
identical  stacks.  Hence,  the  stacking  sequences  of  laminae  in  the  plate  could  be 
defined  by  the  fiber  orientations  of  only  the  top  five  plies.  One  design  variable 
was  assigned  to  each  of  these  five  plies  for  the  representation  of  the  fiber- 
orientation.  The  ply  stacking  sequence  of  the  laminated  plate  was  expressed  as 
[xj,  x,,  X3,  x^,  Xj],/^.  A  layout  of  the  ply  stacking  sequence  of  the  laminated  plate  is 
shown  in  Figure  5.12.  The  design  variables  Xj  through  x^  were  assumed  to  be 
discrete  types,  and  each  was  selected  exclusively  from  the  following  set: 

[-90,  -85,  -80,  -75,     0,  5,  10,     85]  (deg.) 

A  Fourier  series  solution  was  obtained  to  the  mid-point  displacement  of 
the  simply  supported  laminated  plate  under  a  concentrated  load.  A  higher  order 
plate  theory,  which  included  the  effects  of  transverse  shear  deformation  and 
rotatory  inertia  [94]  was  used  in  this  analysis.  Results  were  obtained  for  graphite- 
epoxy  laminated  plates  (E,  =  111. 9  GN/m",  E,  =  10.27  GN/m",  G^,  =  7.31 
GN/m-,  V  =  0.24,  t  =  0.000127  m).  For  the  genetic  search,  a  population  size  of 
50  was  considered.  Each  design  variable,  Xj  through  X5,  was  represented  as  a 
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Figure  5.12  Ply  stacking  sequence  of  laminated  plate. 
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binary  string  of  7  digits.  Four  genetic  searches  were  initially  performed  on  the 
square  plate  with  different  combinations  of  the  probability  of  crossover  P^,  and 
the  probability  of  mutation  P^.  A  population  size  of  50  was  selected  and  the 
genetic  search  was  terminated  after  100  generations.  The  iteration  history  of 
average  central  displacement  for  four  combinations  of  P,  and  P„,  is  shown  in 
Figures  5.13.  The  iteration  history  of  the  smallest  central  displacement  for 
corresponding  value  of  P^  and  P,„  is  shown  in  Figure  5.14.  The  displacement 
shown  in  Figures  5.13  and  5.14  was  nondimensionalized.  This  was  obtained  by 
dividing  the  true  displacement  by  the  displacement  (0.00392m)  of  the  laminated 
plate  with  the  ply  stacking  sequence  as  [-90,  -54,  -18,  18,  54],/,.  The  final  designs 
and  their  displacements  for  these  four  cases  are  shown  in  Table  5.1. 


Table  5.1  Final  designs  and  displacements. 


Genetic  search 
parameters 

Ply-orienlations 

Nondimensionalized 
mid-point  displacement 

Pc  =  0.6,  P„  =  0.01 
P,  =  0.6,  P„  =  0.001 
P,=0.8,  P„,  =  0.01 
Pe=0.8,  P,„  =  0.001 

[-45,  45,-45,  45,-40],/, 
[  40,  -50,  -45,  -45,  45],/, 
[-45,  45,-45,  45,  45],,, 
[-45,  45,  -45,  45,  45],;, 

0.8101 

0.8090 
0.8073 
0.8073 

Another  four  genetic  searches  were  performed  on  the  rectangular  plate 
(aspect  ratio  =  2)  with  the  same  combinations  of  the  probability  of  crossover  P^, 
and  the  probability  of  mutation  P„,.  A  population  size  of  50  was  selected  and  the 
genetic  search  was  terminated  after  100  generations.  The  iteration  history  of 
average  central  displacement  for  four  combinations  of  P,  and  P,„  is  shown  in 
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Figure  5.13  Iteration  history  of  average  displacement  for  the  square  plate. 
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igure  5.14  Iteration  history  of  smallest  displacement  for  the  square  plate. 
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Figures  5.15.  Similar  results  for  the  minimum  central  displacement  are  shown  in 
Figure  5.16.  As  before,  the  displacement  shown  in  Figures  5.15  and  5.16  has  been 
nondimensionalized.  They  were  obtained  by  dividing  the  true  displacement  by  the 
displacement  (0.00183m)  of  a  laminated  plate  with  the  ply  stacking  sequence  as  [- 
90,  -54,  -18,  18,  54]2/3.  The  final  designs  which  had  the  same  values  of  0.8808  for 
displacement  in  these  four  cases  are  shown  in  Table  5.2. 


Table  5.2  Final  designs  and  displacements. 


Genetic  search 
parameters 

Ply-orientations 

Nondimensionalized 
mid-point  displacement 

P,  =  0.6,  P„  =  0.01 
P,=0.6,  P„  =  0.001 
P,=0.8,  P„,  =  0.01 
P,  =  0.8,  P„,  =  0.001 

[  60,  60,  -65,  60,  60]2/, 
[-60,  60,  60,  65,  -60],/, 
[  60,  -60,  -65,  60,  -60]j/, 
[60,  60,  60,  65,  65]y, 

0.8808 
0.8808 
0.8808 
0.8808 

From  the  trend  of  the  numerical  solution  obtained  in  these  two  cases 
(aspect  ratios  equal  to  1  and  2),  the  optimal  ply  orientation  would  be  greater  than 
60  degrees  for  laminated  plates  of  an  aspect  ratio  greater  than  2.  As  an  extreme 
case,  the  optimum  fiber  orientation  would  be  90  degrees  for  all  laminae  as  the 
plate's  aspect  ratio  approaches  infinity.  In  both  these  examples,  it  was  observed 
that  genetic  searches  with  a  higher  crossover  probability,  P,  =  0.8,  located  the  best 
designs  quicker  than  the  search  with  a  lower  crossover  probability  of  P^  =  0.6.  This 
was  attributed  to  the  larger  number  of  new  designs  which  were  created  as  a  result 
of  the  higher  probability  of  crossover.  Genetic  searches  with  the  lower  mutation 
probability,  P„,  =  0.001,  produced  better  average  fitnesses  than  simulations  with  a 
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Figure  5.15  Iteration  history  of  average  displacement  for  the  rectangular  plate. 
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;ure  5.16  Iteration  history  of  smallest  displacement  for  the  rectangular  plate. 
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larger  value  of  Pn,  =  0.01.  This  was  attributed  to  the  lower  degree  of  schema 
disruption  induced  by  the  lower  probability  of  mutation. 

The  real  advantage  to  the  use  of  genetic  search  methods  to  solve  the  ply 
stacking  sequence  problem  becomes  more  evident  when  the  approach  is  compared 
with  proposed  mathematical  programming  solutions  to  the  problem.  Haftka  and 
Walsh  [95]  proposed  a  standout  strategy  for  a  stacking-sequence  design  for  both 
maximum  buckling  load  for  given  total  thickness,  and  minimum  total  thickness 
subject  to  a  buckling  constraint.  The  laminate  was  assumed  to  be  composed 
exclusively  of  0-deg,  90-deg,  and  ±45-deg  plies.  In  order  to  define  the  orientation 
of  each  ply,  four  ply-orientation-identity  variables,  x^,  x^,  X45,  and  x^5  were 
created.  The  variable,  Xq,  x^q,  X45,  or  x^5  is  equal  to  one  if  there  is  a  0-deg,  90-deg, 
45-deg,  or  -45-deg  ply,  respectively.  An  additional  constraint  which  assures  a 
unique  orientation  for  each  ply  is  given  as  follows: 

+  X90  +  X45  +  X_,5  <  1 

A  linear  integer-programming  formulation  of  this  problem  was  solved  using  a 
branch-and-bound  algorithm.  It  was  noted  that  a  large  number  of  sub-problems 
would  be  created  with  the  use  of  such  an  algorithm.  The  number  of  such  sub- 
problems  often  grows  exponentially  with  the  number  of  the  integer/discrete 
variables.  The  computational  expense  would  be  extremely  high  if  a  large  number 
of  possible  fiber  orientations  were  allowed. 
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However,  in  the  genetic  search  method,  only  one  discrete  variable  is 
needed  to  represent  the  fiber  orientation  of  each  ply  regardless  of  the  total 
number  of  allowed  fiber  orientations.  Furthermore,  the  discrete  variables  were 
handled  without  any  more  extra  cost  than  continuous  variables  in  this  approach. 

5.2  Multicriterion  Optimization 
The  multicriterion  optimization  strategies  described  in  preceding  chapter 
were  implemented  for  a  class  of  representative  structural  design  problems.  All 
analysis  was  performed  using  the  finite-element-based  structural  analysis  program 
ANSYS.  The  analysis  was  linked  to  the  genetic  search  procedure  via  pre-  and 
post-processors. 


5.2.1  Statically  Loaded  Ten  Bar  Truss 

A  statically  loaded  ten  bar  truss  shown  in  Figure  5.17,  comprised  the  first 
test  problem.  A  load  P  =  100  kips  was  applied  at  nodes  3  and  5  in  the  direction 
of  displacements  d,  and  dg.  A  truss  bar  length  L  =  360  in.  was  specified.  As 
depicted  in  Table  5.3,  three  variations  of  this  problem  were  constructed,  with 
different  dimensionality  as  represented  by  the  number  of  design  variables  and 
constraints.  The  table  also  shows  the  mix  of  continuous  and  discrete  design 
variables  in  each  problem.  The  sizing  of  this  truss  was  done  with  the  dual  goal  of 
minimizing  both  the  structural  weight  and  the  vertical  displacement  d^.  In 
addition,  displacement  constraints  were  also  imposed  as  specified  in  Table  5.3. 
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Figure  5.17  Ten  bar  truss  structure.  Element  numbers  are  indicated  by  0. 
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Table  5.3  Description  of  the  statically  loaded  ten  bar  truss  problem. 


TEN  BAR  TRUSS  (STATIC  LOAD) 

CASE  1 

Design 
Variable 

Element 

Type 

Constraint 

XI 

1,3 

Continuous 

<  3.0  in. 

X2 

7,8,9,10 

Discrete 

X3 

2,4,5,6 

Discrete 

CASE  2 

Design 
Variable 

Element 

Type 

Constraint 

XI 

1,3 

Continuous 

do  <  3.0  in. 

X2 

2,4 

Discrete 

dj  <  1.5  in. 

X3 

5,6 

Discrete 

X4 

X5 

7,8 

9,10 

Discrete 
Discrete 

CASE  3 

Design 
Variable 

Element 

Type 

Constraint 

XI 

1 

Continuous 

d(,  <  3.0  in. 

X2 

2 

Discrete 

df,  <  1.5  in. 

X3 

3 

Continuous 

dj  <  0.8  in. 

X4-X10 

4-10 

Discrete 

ds  <  0.8  in. 

Elements  of  discrete  set 

0J5  0.50,  0.85,  1.25,  1.90,  2.20,  2.45,  2.60,  3.50,  3.90,  4.60,  4.90,  5.25,  5.75,  6.00  6  55  7  40 
7.70,  8.10,  8.60,  9.00,  9.45,  9.80,  10.25,  10.50,  11.15,  11.65,  12.45,  13.1o,13.70,  14  00  M  45 
14.8M4.90,  15.35,  15.50,  16.50,  17.50,18.75,  20.00  (in^) 
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Different  weights  were  assigned  to  each  of  these  goals.  In  order  to  reduce  the 
differences  between  the  magnitudes  of  individual  criteria,  a  scaling  on  each 
individual  criterion  was  necessary.  The  scaled,  overall  fitness  function  for  this 
problem  was  defined  as  follows,  where  the  scaling  constants  were  determined  by  a 
numerical  assessment  of  the  problem  domain. 


F  =  F 

max 


^3298.9        h.313  * 


(5.11) 


The     term  may  be  thought  of  as  a  penalty  term  introduced  to  account  for  the 
problem  constraints,  much  in  the  manner  of  an  exterior  penalty  function.  is 
the  weighting  factor  for  the  first  criterion,  the  weight  of  the  truss,  and  ranges  from 
0.0  through  1.0  with  a  step  size  0.1.  Similar  variations  were  allowed  for  W„  the 
weight  parameter  associated  with  displacement  d^.  As  defined  in  Table  5.3,  a 
population  size  of  100  was  maintained  at  each  generation,  and  the  genetic 
evolution  stopped  after  40  generations.  Crossover  and  mutation  probabilities  of 
0.8  and  0.01,  respectively,  were  used  in  this  exercise.  Results  were  obtained  for 
each  of  the  defined  test  cases  and  using  both  the  sharing  function  and  the 
vector-evaluated  fitness  strategies. 

Iteration  histories  showing  the  best  fitness  values  corresponding  to  a  given 
weighting  combination  for  cases  1,  2,  and  3  of  Table  5.3,  at  five  different 
generations  of  evolution  are  shown  in  Figures  5.18a-c,  5.19a-c,  and  5.20a-c, 
respectively.  Three,  five  and  ten  design  variables  were  specified  for  cases  1-3.  A 
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gure  5.18a  Best  fitness  values  obtained  using  the  sharing  function  approach. 
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?ure  5.19c  Best  fitness  values  obtained  using  the  VE^^,  approach. 
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comparison  with  results  generated  through  conventional  mathematical 
programming  methods  is  also  shown  on  Figures  5.19a-c.  The  optimal  designs  for 
ten  different  weighting  coefficient  combinations  using  the  sharing  function,  the 
VEi„„  and  the  VE^^^  approaches  for  cases  1,  2,  and  3,  are  summarized  in  Tables 
5.4a-c  through  5.6a-c.  In  general  there  was  good  agreement  among  the  methods 
for  cases  1  and  2  corresponding  to  3  and  5  design  variables;  such  agreement  was 
not  observed  for  the  10  design  variable  problem.  Based  on  our  earlier  discussion 
on  defining  lengths  of^  strings,  it  is  not  difficult  to  account  for  this  trend.  The 
string  lengths  in  the  10  design  variable  case  are  considerably  longer,  allowing  for 
greater  possibilities  of  schema  disruption,  and  hence  also  the  need  to  carry  out 
the  iteration  for  a  larger  number  of  generations.  Approaches  where  such  longer 
strings  may  be  decomposed  into  smaller  defining  length  strings  should  be  explored 
as  possible  remedies  to  the  schema  disruption  problem. 

A  three-criterion  optimization  problem  was  implemented  on  the  statically 
loaded  truss  with  design  variables  and  constraints  defined  in  case  2,  as  shown  in 
Table  5.3.  The  scaled,  overall  fitness  function  for  this  problem  was  a  variation  of 
Eqn.  (5.11)  and  is  defined  as  follows, 


F  =  F 


weight 


(5.12) 


max 


3298.9 


-  1.313 


10,000 


where  the  new  term,  a^.  represents  the  stress  in  the  element  6  of  the  truss.  One 
weighting  variable  which  utilized  a  3-bit  binary  string  was  used  to  form  a  one-to- 
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Table  5.4a  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  1),  using 
sharing  function  approach. 


STATIC-CASE  1  [sharing] 

Weighting 
factors 

Design  variables 

Objective 

Weight 

(ib) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

0.0 

1.0 

19.53 

18.75 

15.35 

0.529 

5844.7 

0.695 

0.1 

0.9 

19.84 

20.00 

10.25 

0.644 

5485.0 

0.697 

0.2 

0.8 

19.84 

18.75 

11.65 

0.760 

5443.4 

0.705 

0.3 

0.7 

19.40 

15.50 

6.55 

0.832 

4320.9 

0.824 

0.4 

0.6 

15.92 

11.65 

5.25 

0.885 

3360.4 

1.045 

0.5 

0.5 

18.21 

9.00 

4.60 

0.887 

2992.2 

1.139 

0.6 

0.4 

14.26 

6.00 

3.50 

0.876 

2163.9 

1.589 

0.7 

0.3 

13.83 

6.55 

4.60 

0.837 

2352.1 

1.478 

0.8 

0.2 

9.95 

4.60 

1.90 

0.711 

1514.7 

2.259 

0.9 

0.1 

6.20 

5.25 

1.25 

0.575 

1332.9 

2.772 

1.0 

0.0 

8.96 

5.75 

2.60 

0.522 

1722.0 

2.003 
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Table  5.4b  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  1),  using  the 
VEjn^  approach. 


STATIC-CASE  1  [VE  J 

Weighting 
factors 

Design  variables 

Objective 

Weight 
(lb) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

0.0 

1.0 

19.95 

20.00 

20.00 

0.501 

6595.0 

0.657 

0.1 

0.9 

19.96 

20.00 

12.45 

0.641 

5740.9 

0.682 

0.2 

0.8 

19.95 

18,75 

8.60 

0.753 

5104.4 

0.727 

0.3 

0.7 

20.00 

14.45 

5.75 

0.828 

4096.2 

0.854 

0.4 

0.6 

19.84 

11.15 

5.25 

0.870 

3502.2 

0.975 

0.5 

0.5 

17.49 

9.00 

5.25 

0.890 

3025.0 

1.134 

0.6 

0.4 

14.79 

8.10 

3.50 

0.871 

2530.0 

1.348 

0.7 

0.3 

11.72 

6.00 

2.60 

0.813 

1918.2 

1.777 

0.8 

0.2 

8.69 

4.60 

1.90 

0.710 

1443.3 

2.364 

0.9 

0.1 

5.97 

4.60 

1.25 

0.557 

1215.8 

2.956 

1.0 

0.0 

5.93 

4.60 

1.25 

0.368 

1213.6 

2.964 
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Table  5.4c  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  1),  using  the 
VE^j.^.  approach. 


STATIC-CASE  1  [VE,J 

Weighting  factors 

Design  variables 

Objective 

Weight 
(lb) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

0.0 

1.0 

19.91 

20.00 

20.00 

0.501 

6592.7 

0.658 

0.1 

0.9 

19.93 

20.00 

13.70 

0.642 

5880.7 

0.676 

0.2 

0.8 

19.89 

18.75 

9.45 

0.754 

5197.2 

0.720 

0.3 

0.7 

19.95 

14.00 

6.55 

0.827 

4111.9 

0.849 

0.4 

0.6 

19.95 

11.65 

5.25 

0.871 

3588.5 

0.953 

0.5 

0.5 

17.86 

9.00 

3.90 

0.887 

2893.2 

1.178 

0.6 

0.4 

15.16 

7.70 

3.50 

0.869 

2486.9 

1.370 

0.7 

0.3 

11.88 

5.75 

2.60 

0.813 

1887.2 

1.806 

0.8 

0.2 

10.05 

4.60 

2.20 

0.711 

1554.3 

2.195 

0.9 

0.1 

6.70 

3.50 

2.20 

0.545 

1188.6 

2.905 

1.0 

0.0 

5.51 

3.90 

2.20 

0.359 

1185.3 

2.986 
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Table  5.5a  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  2),  using  the 
sharing  function  approach. 


STATIC-  CASE  2  [sharing] 

Weighting 
factors 

Design  variables 

Objective 

Weight 
(lb) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

X4 

X5 

0.0 

1.0 

18.SS 

14.85 

6.55 

17.50 

14.00 

0.584 

4801.2 

0.768 

0.1 

0.9 

18.42 

10.25 

4.60 

18.75 

17.50 

0.670 

4784.7 

0.766 

0.2 

0.8 

19.49 

9.45 

4.60 

13.10 

20.0 

0.772 

4547.8 

0.814 

0.3 

0.7 

19.44 

8.10 

5.25 

18.75 

9.45 

0.825 

4113.2 

0.847 

0.4 

0.6 

18.81 

5.75 

2.20 

13.10 

9.00 

0.862 

3285.6 

1.014 

0.5 

0.5 

19.19 

5.25 

4.90 

10.50 

6.55 

0.885 

3025.4 

1.119 

0.6 

0.4 

15.24 

5.25 

4.90 

10.50 

5.75 

0.880 

2737.8 

1.254 

0.7 

0.3 

15.08 

5.75 

2.45 

5.25 

7.70 

0.847 

2354.2 

1.523 

0.8 

0.2 

13.68 

4.90 

4.60 

6.00 

6.00 

0.787 

2272.5 

1.548 

0.9 

0.1 

13.91 

1.90 

0.50 

5.75 

4.60 

0.635 

1751.6 

2.066 

1.0 

0.0 

4.50 

5.25 

0.85 

5.75 

6.55 

0.480 

1584.5 

2.707 
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Table  5.5b  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  2),  using  the 
VEjnj  approach. 


STATIC-CASE  2  [VE,J 

Weighting 
factors 

Design  variables 

Objective 

Weight 
(lb) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

X4 

X5 

0.0 

1.0 

19.96 

20.00 

20.00 

20.00 

20.00 

0.501 

6595.5 

0.657 

0.1 

0.9 

19.96 

16.50 

10.50 

20.00 

20.00 

0.638 

5859.7 

0.672 

0.2 

0.8 

19.93 

11.65 

4.60 

20.00 

13.70 

0.743 

4745.3 

0.747 

0.3 

0.7 

19.67 

9.00 

3.50 

15.35 

11.65 

0.816 

3982.0 

0.851 

0.4 

0.6 

19.67 

9.00 

2.60 

14.00 

9.00 

0.863 

3610.9 

0.931 

0.5 

0.5 

18.78 

5.25 

2.20 

9.80 

8.10 

0.874 

2917.4 

1.135 

0.6 

0.4 

14.46 

4.60 

1.90 

9.00 

6.00 

0.856 

2387.0 

1.386 

0.7 

0.3 

11.04 

4.60 

1.90 

9.00 

5.25 

0.815 

2133.4 

1.584 

0.8 

0.2 

9.16 

1.90 

0.85 

4.60 

4.60 

0.708 

1410.5 

2.405 

0.9 

0.1 

6.36 

1.90 

0.85 

4.60 

4.60 

0.551 

1252.0 

2.750 

1.0 

0.0 

6.43 

1.90 

0.15 

4.60 

4.60 

0.369 

1216.4 

2.990 
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5.5c  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  2),  using  the 
VEg^^  approach. 


STATIC-CASE  2  fVE,J 

Wei 
fac 

ghting 
:lors 

Design  variables 

Objective 

Weight 
(lb) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

X4 

X5 

0.0 

1.0 

19.89 

14.85 

15.35 

20.00 

20.00 

0.512 

6036.9 

0.672 

0.1 

0.9 

19.89 

14.85 

9.00 

20.00 

20.00 

0.638 

5677.5 

0.680 

0.2 

0.8 

19.89 

9.45 

3.90 

20.00 

13.70 

0.743 

4584.0 

0.764 

0.3 

0.7 

19.63 

9.45 

3.50 

15.50 

12.45 

0.817 

4081.2 

0.836 

0.4 

0.6 

19.75 

6.00 

2.60 

13.10 

8.60 

0.859 

3341.6 

0.992 

0.5 

0.5 

17.46 

5.25 

1.90 

10.25 

8.10 

0.874 

2861.7 

1.157 

0.6 

0.4 

14.77 

4.60 

1.90 

8.60 

6.55 

0.856 

2416.6 

1.367 

0.7 

0.3 

12.30 

3.90 

0.85 

6.00 

5.75 

0.807 

1905.6 

1.763 

0.8 

0.2 

9.27 

2.45 

1.25 

5.75 

3.90 

0.700 

1506.5 

2.200 

0.9 

0.1 

7.12 

1.90 

0.85 

3.90 

3.50 

0.534 

1151.0 

2.889 

1.0 

0.0 

7.17 

2.45 

2.20 

3.90 

2.45 

0.357 

1177.3 

2.920 
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Table  5.6a  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  3),  using  the 
sharing  function  approach. 


STATIC-  CASE  3  (shanDgj 

Wc 

fa 

ghcing 

:lors 

Design  vanables 

Obj. 

Weighl 
(lb) 

Disp. 

(m) 

VVI 

W2 

Xl 

.\2 

-\4 

.\f 

X6 

X7 

XS 

.X9 

XIO 

0.0 

1.0 

18.93 

14.45 

18.42 

20.00 

11.15 

11.55 

15J0 

nso 

14.90 

5.75 

0.606 

4824.4 

0.796 

0.1 

0.^ 

1931 

9.00 

16.68 

20.00 

11.65 

4.90 

12.45 

18.75 

14.90 

9.00 

0.700 

4512.8 

0.821 

0.2 

0.8 

1736 

5.75 

18.71 

14.90 

7.70 

2.45 

13.10 

18.75 

11.15 

4.90 

0.785 

38093 

0.909 

OJ 

0.7 

17.13 

2.tO 

19.75 

9.00 

11.65 

3.90 

14.85 

20.00 

12.45 

11.15 

0.846 

4151J 

0.879 

0.4 

O.O 

18.N 

2.45 

18.81 

6.00 

4.90 

3.90 

12.45 

11.15 

1630 

12.45 

0.881 

3641J 

0.961 

OS 

OS 

17.44 

7.40 

14.16 

14.85 

860 

4.60 

12.45 

11.15 

8.60 

4.90 

0.922 

3382.3 

1.076 

0,0 

OA 

W.« 

3.90 

12.07 

13.70 

0.15 

6,55 

13.10 

8.60 

11.15 

4.60 

0.890 

3083.6 

1.082 

0.7 

03 

17.20 

2.60 

18.81 

6.00 

11.65 

030 

4.60 

15  J5 

9.00 

oso 

0.888 

2785.0 

1300 

0.8 

0.2 

5.42 

5.75 

16.83 

1.90 

5.25 

3S0 

10.25 

9.80 

3.90 

4.60 

0.821 

2236.4 

1.827 

0.9 

0.1 

5,29 

5.25 

13.25 

3.50 

4.60 

1.90 

5.75 

6.00 

8.10 

3.90 

0.674 

1906.8 

2.018 

1.0 

0.0 

5.25 

4.60 

5.47 

7.70 

6S5 

0.50 

5.75 

6.00 

8.10 

3.90 

03)6 

18013 

2.408 
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Table  5.6b  Optimal  designs  for  the  statically  loaded  ten  bar  truss  (case  3),  using  the 
VEj,,^  approach. 


STATIC-CASE  3  |VE„^1 

Wc 
fa 

ghling 
•tors 

Design  vanables 

Obj. 

Weighl 
(lb) 

Disp. 
(in) 

Wl 

W2 

Xl 

X2 

.\3 

X4 

X5 

X6 

X7 

X8 

X9 

XIO 

0.0 

1.0 

20.00 

20.10 

19.97 

0.10 

8.10 

20.00 

20.00 

20.00 

2.20 

20.00 

0.659 

4987.4 

0.866 

O.I 

0.9 

20.00 

20.10 

19.97 

0.10 

8.10 

20.00 

20.00 

20.00 

2.20 

20.00 

0.744 

4987.4 

0.866 

0.2 

O.S 

19.W 

19.10 

19.97 

0.10 

810 

20.00 

20.00 

11.05 

1.25 

20.00 

0.826 

4582.4 

0.900 

0.1 

0.7 

19.84 

14.10 

19.97 

0.10 

5.75 

13.70 

20.00 

9.45 

0.15 

20.00 

0.889 

406J.0 

0.975 

O.-l 

0,6 

14.00 

10.10 

19.90 

0.10 

1.25 

9.80 

20.00 

2.45 

030 

1630 

0.922 

3139.6 

1.185 

as 

OS 

12J7 

10.10 

19.98 

0.10 

0.85 

9.80 

18.75 

2.20 

030 

14.45 

0.922 

2948.0 

1,247 

0.6 

0.4 

10.00 

8.10 

HSS 

0.10 

0.85 

8.10 

14.45 

1.90 

0.15 

11.65 

0.901 

2448.2 

1.496 

0.7 

03 

6.89 

8.10 

15.88 

0.10 

0.15 

8.10 

14.00 

0.85 

0.15 

9.45 

0.S40 

20883 

1.738 

O.S 

0.2 

J. 41 

3.10 

15.18 

0.10 

0.15 

3.90 

10.25 

030 

0.15 

9.45 

0.74f. 

1602.3 

2J58 

0.9 

0,1 

5.26 

3.10 

15.03 

0.10 

0.85 

330 

6-55 

030 

0.15 

9.00 

0.602 

1436.2 

2.758 

1.0 

0.0 

6.40 

3.10 

15,03 

0.10 

0.85 

3S0 

4.60 

1.25 

0.15 

9.00 

0.431 

14203 

2.984 
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Table  5.4c  Optimal  designs  for  the  statically  loaded  tenbar  truss  (case  3),  using  the 
VEgj.^  approach. 


STATIC-CASE  3  [VE^,^| 

Wei 
fac 

ghlmg 
lors 

Design  vanables 

Obj. 

(ID) 

Disp. 
(in) 

Wl 

W2 

XI 

X2 

X3 

X4 

X5 

X6 

X7 

X8 

X9 

XIO 

0.0 

1.0 

ISI.W 

20.10 

19.08 

0.10 

15J5 

20.00 

20.00 

20.00 

9.00 

20.00 

0.652 

5427.9 

0.857 

0,1 

0.<J 

19.10 

19.70 

0.10 

14.85 

20.00 

20.00 

20.00 

0.15 

20.00 

0.742 

5048.8 

0.860 

0.2 

0.8 

19A> 

19.10 

19.24 

0.10 

12.45 

18.75 

20.00 

14.00 

0.15 

20.00 

0.830 

4692.4 

0.896 

0.3 

0.7 

l-jJU 

15.10 

18.74 

0.10 

8.60 

15J5 

18.75 

7.70 

0.15 

20.00 

0.898 

4056.3 

o.^»;.' 

0.4 

0.6 

IbAb 

10.10 

18.73 

0.10 

6^5 

11.65 

18.75 

4.60 

0.15 

163 

0.937 

34003 

1.149 

OJ 

O.S 

9.22 

10.10 

19.98 

0.10 

0.15 

9.00 

18.75 

0.85 

0.15 

12.45 

0.921 

2«-2.7 

1J5S 

o.c 

0.4 

9.21 

10.10 

18.73 

0.10 

050 

9.00 

15  J5 

0.50 

0.15 

12.45 

0.897 

2486.8 

1.460 

0.7 

0.3 

5J8 

8.10 

17^ 

0.10 

050 

655 

1535 

0.85 

0.15 

9.45 

0.846 

2111.6 

1.742 

0.8 

0.2 

6.7] 

0.1 0 

10^5 

0.10 

2.20 

6.00 

9.00 

1.90 

0.15 

9.45 

0.751 

1 71 6.4 

2.220 

0.9 

0.1 

2.S4 

5.10 

14.58 

0.10 

0.8S 

7.70 

9.80 

030 

0.15 

4.60 

0.611 

14S7J 

2.r*2 

1.0 

0.0 

2.6S 

6.10 

9J9 

0.10 

2.45 

2.60 

14.00 

030 

0.15 

5.75 

0.44'< 

1482.1 

2.880 
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one  mapping  between  binary  representations  and  the  eight  weighting 
combinations.  The  population  size  at  each  generation  was  50,  and  the  genetic 
search  was  allowed  to  proceed  through  40  generations.  The  iteration  histories  of 
the  best  fitness  value  for  eight  different  combinations  of  weighting  factors  are 
shown  in  Figure  5.21. 

An  interesting  feature  observed  in  the  use  of  the  sharing  function  approach 
was  the  lack  of  a  steady  improvement  in  the  fitness  with  increasing  iterations. 
This  is  largely  attributed  to  the  fact  that  the  binary  bits  representing  the  weighting 
variable  were  altered  during  the  crossover  process,  resulting  in  an  undesirable 
vector-element  switch  and  a  resultant  degradation  in  fitness  value.  An  approach 
where  the  weighting  combination  variable  is  protected  from  the  crossover 
operation  could  be  a  remedy  to  this  problem.  Both  options  in  vector-evaluated 
approach  performed  better  than  the  sharing  approach,  with  the  VE,„^  performing 
slightly  better  than  VE,,^.  Furthermore,  the  only  degradation  encountered  with 
mcreasing  problem  dimensionality  is  the  increase  in  probability  of  schema 
disruption.  Definition  of  shorter  defining  length  strings  is  expected  to  reduce  such 
degradation. 

5-2.2  Dvnamicallv  Loaded  Ten  Bar  Truss 

This  problem  represents  a  variant  of  the  statically  loaded  ten  bar  truss, 
where  a  harmonically  varying  load  f  =  20,000Sin45.2t  N  was  applied  in  the  d, 
direction  at  node  5  as  shown  in  Figure  5.22.  A  truss  bar  length  L  =  100  in.  was 
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W1 
W2 
W3 


Figure  5.21  Best  fitness  values  obtained  using  the  sharing  function  approach  in  the 
3-criterion  problem. 
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20000  sin45 . 2  lb 


Figure  5.22  Dynamically  loaded  ten  bar  truss  structure.  Element  numbers  are 
indicated  by  0. 
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chosen.  A  minimum  weight  design  with  minimum  steady  state  vertical 
displacement  d^  was  sought,  with  constraints  on  dynamic  displacements  and 
member  stresses.  The  design  space  for  this  problem  is  nonconvex  [42];  however, 
it  offers  no  significant  difficulty  to  the  GA  approach.  The  overall  fitness  function 
for  the  problem  was  as  follows: 


F  ^  F 


weight  ,^_^^p 
-  0.493  ^' 


^  98.97 


(5.13) 


As  before,  problems  with  three  different  dimensionalities  were  constructed  (Table 
5.7).  However,  only  results  for  the  sharing  function  approach  and  for  cases  1-3 
are  presented.  Crossover  and  mutation  probabilities  of  0.8  and  0.01,  respectively, 
were  used  in  this  exercise.  W,  and  W,  both  range  between  0.0  through  1.0,  and 
the  condition  Wi  +  W,  =  l  was  enforced.  The  population  size  at  each  generation 
was  100,  and  the  genetic  search  allowed  to  proceed  through  20  generations.  The 
iteration  histories  of  the  best  fitness  value  for  different  weighting  coefficient 
combinations  in  cases  1-3  are  shown  in  Figures  5.23a-c,  and  the  final  designs  of 
the  truss  for  ten  weighting  combinations  for  each  case  are  shown  in  Table  5.8a-c. 

5.2.3  Wing-Box  Structure 

A  wing-box  structure  comprised  the  final  example  with  16  design  variables 
and  38  constraints  [96].  This  problem,  as  shown  in  Figure  5.24,  has  a  combination 
of  stress  and  displacement  constraints;  the  objective  of  the  optimization  problem 
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5.7  Description  of  the  dynamically  loaded  ten  bar  truss  problem. 


TEN  BAR  TRUSS  (DYNAMIC  LOAD) 

CASE  1 

Design  Variable 

Element 

Type 

Constraint 

XI 
X2 
X3 

1,5 
6,7,8,9 
2,3,4,10 

Continuous 
Discrete 
Discrete 

d,  <  0.5  in. 
dg  <  0.5  in. 

CASE  2 

Design  Variable 

Element 

Type 

Constraint 

XI 
X2 
X3 
X4 
X5 

1,3 
2,4 
5,6 
7,8 
9,10 

Continuous 
Discrete 
Discrete 
Discrete 
Discrete 

d,  <  0.5  in. 
dj  <  0.5  in. 
d^  <  0.3  in. 

CASE  3 

Design  Variable 

Element 

Type 

Constraint 

XI 
X2 
X3 
X4-X10 

1 
2 
3 
4-10 

Continuous 
Discrete 

Continuous 
Discrete 

d^  <  0.5  in. 
ds  <  0.5  in. 
dj  <  0.3  in. 

Axial  stresses  (elements  1-10) 
<  45,000  psi 

Elements  of  discrete  set 

0.130,  0.165,  0.175,  0.190,  0.220,  0.245,  0.260,  0.270,  0.280,  0  300  0  315  0  325  0  340  0  3Sn 
0.610,  0.640,  0.670,  0.700,  0.725,  0.765,  0.800,  0.840,  0.880,  0.920,  0.960  1  000  "-''"'^j^;"' 
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(w^-  1-vy,) 

^ure  5.23a  Optimal  designs  for  the  dynamically  loaded  ten  bar  truss  (case  1),  using 
the  sharing  function  approach. 
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w,  (V\^=  1-W^) 

Figure  5.23b  Optimal  designs  for  the  dynamically  loaded  ten  bar  truss  (case  2),  using 
the  sharing  function  approach. 
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(V\^=  1-W^) 

Figure  5.23c  Optimal  designs  for  the  dynamically  loaded  ten  bar  truss  (case  3),  using 
the  sharing  function  approach. 
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Table  5.8a  Optimal  designs  for  the  dynamically  loaded  ten  bar  truss  (case  1),  using 
the  sharing  function  approach. 


DYNAMIC-CASE  1  [sharing] 

Weighting 
factors 

Design  variables 

Objective 

Weight 
(lb) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

0.0 

].0 

0.370 

0.165 

0.840 

0.441 

142.4 

0.217 

0.1 

0.9 

0.300 

0.175 

0.765 

0.440 

131.6 

0.168 

0.2 

0.8 

0.510 

0.175 

0.800 

0.666 

147.4 

0.227 

0.3 

0.7 

0.309 

0.150 

0.400 

0.690 

86.8 

0.301 

0.4 

0.6 

0.210 

0.165 

0.840 

0.772 

133.4 

0.191 

0.5 

0.5 

0.397 

0.165 

0.460 

0.804 

101.0 

0.290 

0.6 

0.4 

0.206 

0.150 

0.390 

0.784 

79.8 

0.370 

0.7 

0.3 

0.102 

0.165 

0.300 

0.766 

66.1 

0.490 

0.8 

0.2 

0.141 

0.165 

0.300 

0.748 

68.4 

0.481 

0.9 

0.1 

0.154 

0.150 

0.340 

0.735 

71.2 

0.430 

1.0 

0.0 

0.102 

0.165 

0.300 

0.668 

66.1 

0.490 
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Table  5.8b  Optimal  designs  for  the  dynamically  loaded  ten  bar  truss  (case  2),  using 
the  sharing  function  approach. 


DYNAMIC-CASE  2  [sharing] 


Wei 
fac 

ghling 
tors 

Design  variables 

Objective 

Weight 
(lb) 

Displacement 
(in) 

Wl 

W2 

XI 

X2 

X3 

X4 

X5 

0.0 

1.0 

0.556 

0.840 

0.670 

0.540 

0.150 

0.003 

172.2 

0.002 

0.1 

0.9 

0.268 

0.340 

0.920 

0.500 

0.340 

0.168 

153.7 

0.007 

0.2 

0.8 

0.204 

0.270 

0.725 

0.315 

0.340 

0.273 

120.3 

0.019 

0.3 

0.7 

0.1  SO 

0.840 

0.175 

0.245 

0.315 

0.398 

112.5 

0.040 

0.4 

0.6 

0.179 

0.765 

0.375 

0.190 

0.150 

0.448 

101.9 

0.030 

0.5 

0.5 

0.141 

0.150 

0.460 

0.435 

0.315 

0.643 

102.5 

0.123 

0.6 

0.4 

0.141 

0.165 

0.920 

0.520 

0.150 

0.763 

123.0 

0.021 

0.7 

0.3 

0.151 

0.350 

0.190 

0.260 

0.500 

0.775 

99.9 

0.113 

0.8 

0.2 

0.159 

0.350 

0.190 

0.150 

0.500 

0.815 

91.6 

0.183 

0.9 

0.1 

0.567 

0.390 

0.340 

0.220 

0.165 

1.028 

104.2 

0.394 

1.0 

0.0 

0.310 

0.375 

0.375 

0.270 

0.165 

0.958 

94.8 

0.263 
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Table  5.8c  Optimal  designs  for  the  dynamically  loaded  ten  bar  truss  (case  3),  using 
the  sharing  function  approach. 


DYNAMIC-CASE  3  [sharing] 

fac 

iors 

Design  variables 

Obj. 

Weight 
(lb) 

Disp. 
(in) 

Wl 

W7 

XI 

X2 

X3 

.\4 

X5 

X6 

X"? 

X8 

X9 

XIO 

0.0 

1.0 

0.129 

0.175 

1.000 

0.450 

1.087 

0.15O 

0.190 

0320 

0.960 

0.640 

0.499 

171.4 

0.246 

0.1 

0.9 

0J20 

0.8«0 

0.610 

03Z5 

0.974 

0.175 

0.175 

0.800 

0.165 

0320 

0.668 

161.0 

0.277 

0.2 

0.8 

0.753 

0315 

0.435 

0.245 

1.070 

0.150 

0315 

0.670 

0380 

0350 

0.809 

1583 

0301 

03 

0.7 

0.749 

0.670 

0.270 

0.700 

0J57 

0.150 

0325 

0.435 

0360 

0.400 

0.905 

1533 

0..309 

0.4 

0.6 

1.020 

0390 

0.400 

0.460 

0305 

0.150 

0350 

0.450 

0300 

0.190 

0.907 

1283 

0319 

O-'i 

0-S 

0.522 

0.420 

0.S40 

0380 

0392 

0315 

0.270 

0.190 

0.175 

0315 

1.079 

124.9 

0.442 

0.6 

0.4 

0.439 

0.485 

0340 

0.270 

0J29 

0.150 

0340 

0315 

0.640 

0.15O 

1.034 

1203 

0375 

0.7 

03 

0341 

0.400 

0300 

0.245 

0.289 

0.175 

0375 

0.245 

0.485 

0.165 

0.975 

1003 

0.434 

O.S 

0.2 

0.424 

0375 

0^80 

0.460 

0313 

0.175 

0315 

0.340 

0.420 

0350 

1.180 

1263 

0387 

0.9 

0.1 

0.936 

0.150 

0350 

03O0 

0391 

0.150 

0.450 

0340 

0.670 

0.150 

1.29.1 

KVl3 

0.353 

1.0 

0.0 

0374 

0315 

0.220 

0.2*0 

1.106 

0.190 

0300 

0.450 

0340 

0.640 

1.438 

1423 

0.435 
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Figure  5.24  Statically  loaded  wing-box  structure. 
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was  to  achieve  the  minimum  weight  of  the  wing-box  structure  and  a  maximum 
value  of  the  sum  of  the  first  two  natural  frequencies.  Problem  parameters  used  in 
the  numerical  study  are  described  in  Table  5.9.  The  overall  fitness  function  was 
defined  as  follows: 


weight  ^ 


500 


0.5xJ2_.0.5x.l50 


freq. 


(5.14) 


The  population  size  was  maintained  at  100  and  the  search  was  carried  out  to  30 
generations.  Sample  results  for  both  the  sharing  function  and  vector-evaluated 
approaches  are  presented  here.  The  evolution  of  the  best  fitness  for  different 
weighting  combinations  are  shown  in  Figures  5.25a-c.  These  trends  are  similar  to 
those  obser\'ed  in  the  smaller  dimensionality  problems.  A  comparison  with  results 
obtained  using  nonlinear  programming  methods  show  that  a  larger  number  of 
evolution  cycles  are  required  for  convergence  in  both  the  sharing  function  and  the 
VEg^.^  approaches.  The  results  are  very  close  to  converged  values  when  using  the 
^Einc  approach.  In  general,  the  current  implementation  of  the  sharing  function 
approach  was  consistently  outperformed  by  the  vector-evaluated  approach. 
Additionally,  there  is  a  trend  that  the  VE-^^  tends  to  outperform  the  VE^^, 
approach,  particularly  with  increasing  problem  dimensionality. 

It  is  difficult  to  compare  the  computational  efficiencies  of  the  gradient- 
based  nonlinear  programming  method  and  the  genetic  search  method.  A 
nonlinear  programming  algorithm  will  iterate  until  a  relative  minimum  or 
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Table  5.9  Description  of  the  statically  loaded  wing-box  structure. 


WING-BOX 

Design 

V  ul  IclUlC 

Element 

Element 
Type 

Node  Connectivity 

Constraint 

VI 

X2 

1,2 
3,4 

1 E 
TE 

1-3,  8-10 
3-5,  10-12 

U,  <  2.0  in. 
for  all  unconstrained  nodes 

X3 

5,6 

TE 

2-4,  9-11 

X4 
X5 

7,8 

9,10 

TE 
TE 

4-6,  11-13 
6-7,  13-14 

Axial  or  equivalent 
stresses  <  10,000  psi 
for  elements  (1-28) 

X6 

11,12,13,14 

TM 

1-2-4,  4-3-1, 
8-9-11,  11-10-8 

X7 

15,16,17,18 

TM 

3-4-6,  6-5-3, 
10-11-13,  13-12-10 

X8 

19,20 

TM 

5-6-7,  12-13-14 

X9 

21 

QM 

1-8-10-3 

XIO 

22 

QM 

3-10-12-5 

Xll 

23 

QM 

2-9-11-4 

TE:  Truss  element 

X12 

24 

QM 

4-11-13-6 

TM:  Triangular  membrane 

X13 

25 

QM 

6-13-14-7 

QM:  Quadrilateral 
membrane 

X14 

26 

QM 

3-10-11-4 

X15 

27 

QM 

5-12-13-6 

0.01  <  X,.,  <  3.50  in' 

X16 

28 

QM 

5-12-14-7 

0.01  <  X^,^  <  3.50  in 
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■      l8t  oeneration 

+  10th  generation 

o  20th  o«r>eration 

A  30th  generation 

V  Nonlrear  pr oyaming 


Figure  5.25b  Best  fitness  values  obtained  using  the  VEj^,  approach. 
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Figure  5.25c  Best  fitness  values  obtained  using  the  VE^^^  approach. 


162 

maximum  is  reached,  regardless  of  how  many  function  evaluations  are  necessary. 
Genetic  search  algorithms  either  finish  when  convergence  to  a  best  solution  is 
reached,  or  when  the  number  of  fitness  function  evaluations  reaches  a  prescribed 
maximum  cap.  In  genetic  search,  the  optimum  is  usually  obtained  before  the 
entire  population  converges  to  the  optimal  design.  It  is,  therefore,  practical  to 
perform  genetic  searches  by  assigning  a  maximum  value  to  the  number  of  fitness 
function  evaluations  which  are  computationally  affordable. 

An  interesting,  though  indirect  comparison  can  be  made  between  the 
efficiencies  of  these  two  approaches,  which  have  totally  different  convergence 
criteria.  For  the  five-variable  case  of  the  statically  loaded  ten-bar  truss  problem, 
it  took  an  initial  total  of  2,748  function  evaluations  to  perform  11  optimizations, 
based  on  different  weighting  combinations.  More  function  evaluations  will  be 
required  when  the  optimization  problem  was  re-run  with  different  starting  design 
estimates.  This  procedure  is  recommended  in  working  with  gradient  based 
methods,  as  the  existence  of  relative  minima  in  the  design  space  is  not  generally 
known  apriori.  Eleven  'globally  optimal'  solutions  obtained  are  shown  in  Figures 
5.19a-5.19c.  As  a  comparison,  4,000  fitness  function  evaluations  were  necessary  in 
the  genetic  search,  and  the  results  obtained  by  the  VE,„,  and  VE,,,  approaches 
were  observed  to  be  as  good  as  these  'globally  optimal'  designs.  It  is  noted  that 
only  four  discrete  design  variables  are  involved  in  this  optimization  problem. 
When  a  large  number  of  discrete/integer  design  variables  are  presem,  the  branch- 
and-bound  approach  used  in  nonlinear  programming  techniques  to  accoum  for  the 
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non-continuous  design  variables  requires  many  more  function  evaluations.  Genetic 
search  algorithms  have  no  additional  computational  requirements  due  to  an 
increased  number  of  discrete/integer  design  variables  in  the  optimization 
problem. 


5-2.4  Viscoelasticallv  Damped  Laminated  Composite  Re,?^m 

An  integrated  approach  to  the  multicriterion  optimal  design  of  a 
viscoelastically  damped  aluminum  cantilever  beam  was  investigated  in  Chapter  3. 
To  investigate  the  discontinuous  application  of  damping  material  on  the  beam,  a 
genetic  search  approach  was  used.  The  results  described  below  were  obtained  by 
considering  a  population  size  of  100.  A  combination  of  the  VE^^,  and  VE„,^ 
strategies  was  used  to  account  for  the  multicriterion  nature  of  the  problem  and 
values  of  P,  =  0.8  and  P„,  =  0.01  were  used.  In  this  problem,  the  base  structure  was 
a  graphite-epoxy  composite  laminate  with  a  layup  sequence  given  as  [0,  45,  90,  - 
45,  0]y^.  The  material  properties,  as  shown  in  Table  5.10,  were  used  in  this 
simulation.  The  beam  was  divided  into  20  equal  segments  along  its  length,  and  a 
0/1  type  of  variable  was  used  to  denote  the  absence/presence  of  damping 
material  on  a  given  segment.  This  resulted  in  a  20-digit  string  of  O's  and  I's 
defining  the  placement  of  the  damping  material.  The  structural  system  of  the 
composite  sandwich  beam  is  shown  in  Figure  5.26.  The  design  variable  set 
included  variables  indicating  the  presence  of  damping  material  on  a  beam 
segment  X^^^  [Xj,  X,,     X20],  the  thickness  of  the  viscoelastic  layer  X^^,  the 
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Table  5.10  Material  properties  of  the  structural  system. 


Constraining 
layer 

Viscoelastic 
layer 

Graphite-epox}' 
lamina 

(GN/m-) 

69.0 

0.0021 

127.9 

(GN/m-) 

10.27 

0.3 

0.499 

0.24 

0.20 

Damping 
factor 

0.0 

0.10 

0.005 

Density 
(kg/m^) 

2700.0 

970.0 

1580.0 
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Figure  5.26  Viscoelastically  damped  laminated  composite  beam  structure. 
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thickness  of  the  constraining  layer  X^^n'       the  ply  thickness  of  the  base  structure 
Xpiy.  Here,  X^^^,  X^^  and  Xp|y  were  discrete  variables,  selected  from  lists  of 
allowable  values  shown  in  Table  5.11.  The  finite  element  analysis  described  in 
Chapter  3  was  used  to  evaluate  the  structural  damping  and  natural  frequencies,  as 
well  as  the  static  displacements  that  resulted  when  different  types  of  damping 
treatments  were  applied. 

The  vector  of  criterion  contained  four  elements;  the  first  three 
corresponded  to  the  damping  ratio  in  the  first  three  modes  and  the  last  was  the 
total  weight  of  the  beam.  The  constraint  set  comprised  lower  and  upper  bounds 
on  the  first  three  natural  frequencies  and  a  maximum  allowable  displacement 
under  a  static  load  of  1-N  applied  at  the  tip  of  the  beam.  These  design 
constraints  are  summarized  in  Table  5.11.  The  mathematical  statement  of  the 
optimum  design  problem  can  be  written  as  follows: 


Minimize 


damp^       -  damp^_       ^  damp^  ^weight"'' 


(5.15) 


Subject  to        gj  <  0,  y=l,2,...,7 


where  damp^  represent  the  i'^mode  damping,  and  weight  represents  the  total 
weight  of  the  structural  system.  The  superscript  av  represents  the  populational 
average  value  of  the  system  variable  in  a  given  generation.  A  slowed  searching 
rate  in  the  genetic  algorithm  may  result  from  the  presence  of  significantly 
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Table  5.11  Design  constraints  and  allowable  discrete  values  of  design  variables 
^ly  Xcon,  and  X^,. 


.000150  .000175  .000200  .000225  (m) 

V 

con 

.0001  .0002  .0003  .0004  .0005  .0006 
.0007  .0008  .0009  .0010  .0011  .0012 
.0013  .0014  .0015  .0016  .0017  .0018 
.0019  .0020  .0021  .0022  (m) 

.0001  .0002  .0003  .0004  .0005  .0006 
.0007  .0008  .0009  .0010  (m) 

Design 
Constraints 

20.0  <  freqi  <  50.0  (Hz) 
130.0  <  freq,  <  200.0  (Hz) 
400.0  <  freq3  <  700.0  (Hz) 
tip  vertical  displacement  <  0.002  m 
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damping  values  in  certain  designs.  These  low  values  would  most  probably  be 
present  only  in  early  generations.  In  order  to  avoid  such  situations,  the  average 
value  of  a  population,  instead  of  a  constant  value,  was  used  to  scale  the 
magnitudes  of  these  criteria.  Wj,  1  =  1,4,  represent  weighting  parameters  and  each 
weight  was  used  to  emphasize  the  relative  importance  of  the  i"*  criterion.  To 
account  for  the  problem  constraints,  a  penalty  term     was  introduced  as  an 
exterior  penalty  function.  The  scaled,  overall  fitness  function  for  this  probl  em  was 
defined  as  follows: 


F  =  F 

max 


damp^       -  damp 2       ^  damp^        weight'^'  * 


(5.16) 


For  five  different  weighting  combinations  of  these  candidate  criteria,  the 
optimal  designs  obtained  from  the  genetic  search  procedure  are  presented  in 
Table  5.12.  It  is  noted  replacing  the  aluminum  beam  with  a  graphite-epoxy 
laminated  beam  gave  a  50%  increase  in  damping  and  curtailed  the  weight  by  30% 
m  this  application.  According  to  the  results,  a  continuous  application  of  damping 
tapes,  which  spans  from  the  root  of  the  beam  to  a  point  close  to  the  end,  is  a 
reliable  selection  for  obtaining  high  damping  values  in  all  three  modes.  For 
maximum  damping  in  the  2"'^  and  3^'^  modes,  the  most  effective  region  of  damping 
tape  application  was  near  the  middle  area.  Final  designs  with  discrete  damping 
treatments  were  obtained  for  the  weighting  combination  in  which  the  first-mode 
damping  was  weighted  with  0.5,  and  for  the  weighting  combination  wherein  the 
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Table  5.12  Multicriterion  optimal  designs. 


Weighting 

combination 

uam  \     1      /u/  L  oJ 

Loss  Factor 

W  T 

W, 

Jply 

Weight 

V 

con 

Mode  1 

(N) 

W2 

Mode  2 

W3 

(m) 

Mode  3 

W 

0.5 

11111111111101100000 

0.0230 

U.U 

.000200 

0.0155 

1.2415 

0.0 

.0022 

0.0149 

0  5 

.UUUl 

U.U 

11011111111111111000 

0.0135 

U.J 

.000200 

0.0266 

1.3125 

0.0 

.0022 

0.0197 

0.5 

.0001 

0.0 

00111111111111111110 

0.0136 

0.0 

.000200 

0.0250 

1.3481 

0.5 

.0022 

0.0233 

0.5 

.0001 

0.2 

11111111111111111100 

0.0221 

0.2 

.000200 

0.0252 

1.3836 

0.2 

.0022 

0.0230 

0.4 

.0001 

0.1 

11111111111111111000 

0.0196 

0.1 

.000200 

0.0216 

1.1860 

0.1 

.0016 

0.0185 

0.7 

.0001 
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second-mode  was  weighted  with  0.5.  The  non-continuous  damping  treatment 
configurations  used  in  these  two  designs  would  not  be  feasible  for  optimization 
using  the  two  continuous-type  models  presented  in  Chapter  3.  Gradient-based 
nonlinear  programming  methods  are  essentially  inept  in  solving  the  problems 
presented  in  this  section  due  to  the  presence  of  the  large  numbers  of  discrete 
design  variables. 


CHAPTER  6 


APPROXIMATION  CONCEPTS  IN  GENETIC  ALGORITHM  BASED 

OPTIMIZATION 


6.1  Introduction 

A  definitive  gap  can  be  identified  between  a  generally  efficient,  gradient- 
based  mathematical  programming  algorithm,  and  a  computationally  intensive 
direct  enumerative  scheme.  As  discussed  earlier,  the  former  approach  is  the 
method  of  choice  in  problems  where  the  design  space  is  continuous  and  convex. 
It  IS,  however,  prone  to  converging  to  a  relative  optimum  when  the  design  space  is 
nonconvex.  While  not  adversely  affected  by  nonconvexities,  enumerative  schemes 
place  a  severe  requirement  on  the  computational  resources.  Genetic  algorithms 
represent  a  relatively  new  approach  for  bridging  this  gap.  They  inherently  possess 
a  better  chance  of  locating  the  global  optimum  than  nonlinear  programming 
algorithms,  and  though  computationally  intensive,  require  less  effort  than 
enumerative  search  techniques.  New  genetic  search  strategies  having  significantly 
improved  efficiency  may  be  available  in  the  future;  however,  present  methods 
usually  require  thousands  of  fitness  function  evaluations  to  achieve  an  optimum 
solution  for  a  structural  optimization  problem  having  only  moderate 
dimensionality.  More  computational  effort  is  needed  when  the  optimizati 
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problem  involves  a  large  and  complicated  design  space.  Computationally 
expensive  finite  element  analysis  methods  are  often  used  to  evaluate  structural 
response  during  an  optimal  design  cycle.  Since  genetic  algorithms  require  that  this 
analysis  be  performed  several  times,  it  is  often  prudent  to  replace  the 
computationally  intensive  analysis  by  a  less  demanding  approximation.  Hence, 
approximation  concepts  in  genetic  algorithms  are  a  subject  worthy  of  investigation. 
A  computationally  efficient  alternative  to  the  finite  element  analysis  method, 
which  also  yields  an  acceptable  level  of  accuracy,  is  desirable  for  use  in  genetic 
algorithms.  The  applicability  of  using  an  approximation,  based  on  a  trained 
neural  network,  is  proposed  as  an  alternative  to  the  finite  element  analysis 
method.  The  problem  under  consideration  is  one  pertaining  to  the  design  of  a 
viscoelastically  damped  composite  beam  structure. 

6.2  Neural  Networks 
The  neural  computing  concept  has  its  philosophical  basis  in  the  biological 
cerebral  process.  Neural  networks  are  parallel  computing  models  which  use  many 
small  elements,  referred  to  as  neurons  or  nodes,  that  are  densely  interconnected. 
A  set  of  inputs  to  the  network  relays  through  the  neurons  in  a  pattern  that  can 
alter  network  characteristics  in  a  learning  phase,  or  be  subsequently  recognized  in 
a  memory  operation.  It  has  been  shown  [97]  that  a  neural  network  can  model  the 
behavior  of  any  linear  or  nonlinear  continuous  system.  The  use  of  neural  networks 
in  the  analysis  of  structural  systems  holds  promise  as  a  new  tool  for  efficiently 
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solving  problems  that  may  be  computationally  intensive  or  hold  solutions  that  are 
not  readily  apparent  by  traditional  methods.  Recent  studies  in  this  area  have  been 
performed  by  Hajela  [98],  where  a  network  was  trained  to  perform  automated 
structural  design.  A  similar  approach  was  used  by  Rehak  [99],  who  trained  a 
network  to  represent  the  stiffness  matrix  for  a  two  degree-of-freedom  spring 
system.  Neural  networks  have  been  trained  to  model  the  force-displacement 
relationship  in  a  number  of  simple  structural  systems.  The  applicability  of  using 
neural  networks  in  genetic  search  optimization  problems  to  model  complex 
relationships  between  input  and  output  of  the  structural  system,  is  examined  in 
this  work. 

The  configuration  of  a  neural  network  is  defined  in  terms  of  individual 
neurons,  the  network  connectivity,  the  weights  associated  with  interconnections 
between  neurons,  and  the  activation  function  for  each  neuron.  The  network  maps 
a  set  of  input  vectors  to  a  set  of  output  vectors.  With  appropriate  training,  the 
neural  network  can  "learn"  to  perform  reasonably  accurate  mapping.  The 
functional  operation  of  the  neuron  is  first  described. 

Consider  the  single  neuron  shown  in  Figure  6.1.  This  neuron,  referred  to  as 
neuron  q,  receives  a  set  of  m  inputs,  x^,  j  =  l,m  and  produces  an  output  N^.  There 
is  a  weight  associated  with  each  input.  The  input-output  relation  is  obtained  as 


N  =  F{A)  =  F 


f 

m 


(6.1) 
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Figure  6.1 


A  single  neuron. 
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where  A  is  the  weighted  sum  of  all  inputs  comiected  to  the  neuron  q,  and  the 
function,  F,  is  the  activation  function  which  transforms  the  weighted  sum  on 
inputs  to  an  output  response.  This  activation  function  can  be  the  source  of 
nonlinearity  in  the  neuron,  and  in  the  present  work,  a  sigmoid  function  was  used 
for  this  purpose.  The  sigmoid  function  is  defined  as, 

=  (6.2) 

• ,  r 

where  B  is  the  bias  parameter  used  to  modulate  the  output  of  the  neuron.  The 
principal  advantage  of  using  the  sigmoid  function  is  its  ability  to  handle  both  large 
and  small  numbers.  Referring  to  Figure  6.2,  it  can  be  seen  that  the  outputs  for 
very  large  and  very  small  numbers  approach  1  and  0,  respectively,  while  numbers 
in  between  have  smooth  variations.  The  output      of  the  neuron,  obtained  from 
the  sigmoid  function  is  then  treated  as  an  input  for  neurons  in  next  layer. 

The  neural  network  of  this  study  used  an  input  layer  consisting  of  a 
number  of  neurons  that  receive  design  variables  values  of  the  structural  systems, 
and  an  output  layer  consisting  of  a  number  of  neurons  that  produce  values 
representing  the  resultant  behavioral  response  of  the  structural  system.  If  no 
hidden  layers  exist  between  the  input  and  output  layers,  the  network  is  referred  to 
as  a  flat  network.  Simple  flat  networks  show  limited  applicability  in  the 
representation  of  nonlinear  relationships.  A  neural  network  with  a  number  of 
hidden  layers  is  referred  to  as  hidden-layer  network.  Each  hidden  layer  may  have 


176 


STiall  (A-B)  0  large  (A-B) 


Figure  6.2  A  typical  sigmoid  function. 


177 

an  arbitrary  number  of  neurons.  Sketches  of  a  flat  network  and  two  different 
hidden-layer  networks  are  shown  in  Figure  6.3.  A  flat  network  which  has  a  3- 
neuron  input  layer  and  a  2-neuron  output  layer  is  referred  to  as  having  a  (3-2) 
architecture.  The  same  network  with  a  2-neuron  hidden  layer  is  described  as 
having  (3-2-2)  architecture;  likewise,  the  presence  of  two  2-neuron  hidden  layers 
would  be  denoted  by  a  (3-2-2-2)  architecture.  After  a  neural  network  is  built,  it 
must  be  trained  to  produce  an  acceptable  output  vector  for  a  given  set  of  input 
vectors.  One  supervised  learning  technique  used  for  this  purpose  is  the  back- 
propagation  training  algorithm  [100].  This  is  especially  well  suited  for  neural 
networks  with  multilayer  architectures. 


6.3  Back-Propagation  Training  Algorithm 
Prior  to  training  a  neural  network,  the  input  and  output  vectors 
representing  the  structural  system  parameter  values  must  be  determined,  either 
analytically  or  by  numerical  methods.  Randomly  generated  small  values  are 
initially  assigned  to  all  network  weights  and  bias  parameters.  The  first  output 
vector  of  the  network  is  computed  from  the  input  vector  using  these  initial  weights 
and  the  bias  parameters.  The  error  between  an  element  of  the  output  vector  and 
the  expected  output  at  that  element  is  computed  as 

=      -  (6.3) 


where     is  the  network  output  of  the  neuron  i,  as  defined  in  (6.1),  and  T,  is  the 
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re  6.3  Representative  neural  network  architectures. 
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expected  output  of  the  same  neuron.  This  error  can  be  reduced  by  changing  the 
values  of  weights  connecting  all  neurons  in  the  previous  layer  to  the  neuron  that  is 
under  examination.  Weight  adjustments  are  computed  from  the  following 
equation:  , 

^^P„  =  y^.,N^j  '    '     '  (6-4) 

where  Npj  denotes  the  output  of  neuron,  p,  in  the  hidden  layer,  j,  which  directly 
precedes  the  output  layer,  k;  AWp,^  is  the  change  in  value  of  the  weight  between 
the  neuron,  p,  in  the  previous  hidden  layer  and  the  neuron,  i,  in  the  output  layer 
k;  Y  represents  the  learning  rate  coefficient.  The  learning  rate  coefficient,  typically 
selected  between  0.2  and  0.9,  functions  as  a  step  size  parameter  in  the  back- 
propagation  process,  and  determines  the  speed  of  the  network  training;  6,^  is  the 
product  of  the  error  signal  defined  in  (6.3)  and  the  derivative  of  the  output  N  of 
the  sigmoid  function  F(A)  of  the  present  neuron,  and  is  obtained  as  follows: 

''^•  =  -ar^'>=^'^(^  ^^-^^ 

where  the  subscripts  i  and  k  denote  the  i"^  neuron  in  k"^  layer  (k"'  layer  is  the 
output  layer  of  the  reural  network  here).  A  modification  of  (6.4)  is  proposed  by 
Rumelhart  et  al  [101],  who  added  a  momentum  term  into  (6.4)  to  stabilize  the 
training  and  incorporate  the  memory  into  the  adjustment  of  the  weights.  This   '  ,  " 
modifies  (6.4)  as  follows. 
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where  superscript  t  represents  the  cycle  of  the  iteration,  and  a  represents  the 
momentum  coefficient. 

The  error  signal,  E„  is  obtainable  in  the  form  defined  in  (6.3)  only  for 
neurons  in  the  output  layer  because  expected  output  T,  is  unknown  for  the 
previous  layers.  For  neurons  in  the  hidden  layer  adjacent  to  the  output  layer,  the 
following  equation  is  used  to  replace  (6.5). 

-         E  ^M- (6.7) 

i 

where  6^^  is  the  <5  corresponding  to  the  p"^  neuron  in  the  hidden  layer  j.  After  8 
is  obtained  for  neurons  in  this  hidden  layer,  the  adjustment  of  weights  between 
neurons  in  the  next  hidden  layer  and  the  neuron  under  consideration  is  obtained 
using  (6.4)  or  (6.6)  with  appropriate  changes  in  indices.  This  process  is  repeated 
for  all  hidden  layers.  The  back-propagation  training  process  repeatedly  considers 
all  training  patterns  umil  the  error  is  reduced  to  less  than  a  prescribed  value. 
Another  approach  involves  the  consideration  of  the  accumulated  squared  error 
sum  from  all  the  training  input-output  vectors.  The  training  process  is  then  itself 
an  error  minimization  problem. 
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6-4  Miilticriterion  Optima]  Design  of  Three-Dimensional  Viscoelastically  Damped 

Laminated  Composite  Beam 

A  genetic  search  approach  was  implemented  for  the  optimal  design  of  a 

three-dimensional,  viscoelastically  damped,  laminated  composite,  cantilevered 

tapered  beam.   The  design  required  a  simultaneous  minimization  of  the 

structural  weight  and  maximization  of  structural  damping.  A  diagram  detailing 

the  structure  of  the  beam  is  shown  in  Figure  6.4.  The  structural  system  comprised 

of  a  20-ply,  unidirectional  graphite-epoxy  laminated  composite  base  structure  on 

which  viscoelastic  damping  treatment  was  applied.       Each  ply  had  a  uniform 

thickness  that  was  determined  by  a  discrete  design  variable  Xjg.  This  design 

variable  assumes  a  value  selected  from  the  set  of  [0.000100,  0.000125,  0.000160, 

0.000225]  (m).  The  composite  laminate  was  symmetric  and  balanced;  therefore, 

the  stacking  sequence  could  be  represented  by  the  five  design  variables,  [x^.,  x,^, 

^i5»  Xjo,  ^vh/s-  These  five  variables  were  defined  as  discrete  variables  that 

represented  the  fiber  orientations  of  the  top  five  plies  of  the  laminate.  Each  of 

the  five  design  variables  was  selected  from  the  following  set  [-90,  -85,  -80,  80, 

'  85]  (degrees).  Design  variables  which  defined  the  composite  base  beam  and  the 

stacking  sequence  of  20  graphite-epoxy  laminae  are  shown  in  Figure  6.5. 

The  top  surface  of  the  composite  beam  was  partitioned  into  twelve 

sections.  Six  sections  of  the  same  length  were  distributed  on  each  side  of  the 

cantilevered  tapered  beam.  The  damping  treatment  of  the  beam  was  specified  by 

twelve  design  variables,  x^  through  x,^,  of  0/1  type.  The  damping  treatment 

method  is  illustrated  in  Figure  6.4.  Each  of  the  12  sections  could  have  the 
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Figure  6.4  Viscoelastically  damped  laminated  composite  beam. 


Figure  6.5  Ply  sequence  arrangement  of  laminated  beam. 
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damping  tape  independently  applied.  When  two  consecutive  sections  were  both 
subjected  to  damping  treatment,  a  one-piece  damping  tape  was  applied.  The  ply 
thickness  of  each  graphite-epoxy  laminae,  Xjg,  the  thickness  of  the  constraining 
layer,  Xjg,  and  the  thickness  of  the  viscoelastic  layer,  Xjq,  were  all  defined  as 
discrete  quantities  with  admissible  values  selected  from  the  set  shown  in  Table 
6.1(a). 

As  stated  earlier,  the  dual  goal  of  this  optimization  problem  was  to  both 
maximize  the  damping  of  the  first  three  vibrational  modes  and  minimize  the 
structural  weight.  The  design  constraint  set  comprised  of  the  lower  bounds  on  the 
first  three  natural  frequencies,  and  the  maximum  displacement  when  a  static  load 
of  1-N  was  applied  at  the  tip  of  the  beam.  The  maximum  work  theory  [102], 
referred  to  as  the  Tsai-Hill  theory,  was  used  as  a  failure  criterion  on  the  graphite- 
epoxy  composite  laminate.  The  external  loading  condition  used  in  the  failure 
criterion  was  the  combination  of  N,=250,000  N/m  and  M,  =  150  Nm/m.  Material 
properties  of  the  structural  system  and  design  constraint  specifications  are  shown 
in  Table  6.1(a)  and  (b).  The  scaled,  overall  fitness  function  was  expressed 


as 


F  =  F 

max 


damp^         damp,        damp^  ^,p;„Uf 
[     damp^      -  damp,        damp^      '  weight ' 


(6.cS) 


where  w„  i-1,4,  represents  the  weighting  parameters  used  to  emphasize  the 
relative  importance  of  the  i"'  criterion.  For  the  calculation  of  the  penalty  term,  P„ 
the  following  scaled  quantities  need  to  be  defined. 
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Table  6.1(a)  Constants  and  discrete  sets  of  design  variables  Xjg,  Xjg,  and  X. 


^18 

.000100  .000125  .000160  .000225  (m) 

Xjg 

.0001  .0002  .0003  .0004  .0005 
.0006  .0007  .0008  .0009  .0010  (m) 

.0001  .0002  .0003  .0004  .0005 
.0006  .0007  .0008  .0009  .0010  (m) 

Design 
constraints 

1)  4.0  <  freqi  (Hz) 

2)  20.0  <  freq,  (Hz) 

3)  50.0  <  freqj  (Hz) 

4)  tip  vertical  displacement  <  0.05  m 

5)  failure  criterion  with  1 

=  250kN/m,  and      =  150Nm /m  | 

Table  6.1(b)  Material  properties  of  the  structural  system 


constraining 
layer 

Viscoelastic 
layer 

Graphite-epox7 
lamina 

El 

(GN/m-) 

69.0 

0.0021 

127.9 

E, 

(GN/m-) 

10.27 

0.3 

0.499 

0.24 

^2 

0.20 

Damping 
factor 

0.0 

0.10 

0.005 

Density 
(kg/m^) 

2800.0 

970.0 

1580.0 
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5.0 
freq^ 

20.0  , 
freq2 

s.-^-l  (6.9) 

^  005  ' 

^5  =  (  — )    -  (  —  )(  —  )  -  {-^)    -  (^)     -  1 


where  g^,  i  =  1,5,  represents  the  scaled  value  of  the  i""  constraint.  A  design  is 
feasible  when  these  quantities  are  all  negative.  These  constraint  values  were  used 
directly  in  the  input-output  mapping  available  in  the  neural  network  approach.  A 
total  of  twenty  design  variables  and  eight  output  quantities  were  involved  in  the 
optimization  problem.  These  eight  output  quantities  consisted  of  damping  in  the 
first  three  modes,  and  five  previously  defined  scaled  constraint  quantities,  g,.  An 
alternative  to  this  arrangement  would  be  to  use  actual  values  of  natural 
frequencies  and  displacement  as  elements  of  the  output  vector  of  the  network. 

In  order  to  explore  the  effectiveness  of  different  architectures  in  this 
problem,  a  number  of  neural  network  architectures  were  implemented.  An  initial 
model  consisted  of  an  input  layer  of  20  neurons,  two  hidden  layers,  each 
consisting  of  30  neurons,  and  an  output  layer  with  8  neurons.  The  twenty  neurons 
in  the  input  layer  represented  the  twenty  design  variables  of  the  structural  system. 
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The  eight  neurons  in  the  output  layer  represent  the  eight  outputs  of  the  structural 
system  discussed  earlier.  This  neural  model  is  referred  to  as  (20-30-30-8),  and 
was  used  to  represent  the  input-output  relationship  of  the  structure.  A  number  of 
smaller  models  were  also  used  as  an  alternative  to  the  (20-30-30-8)  architecture 
for  the  present  problem.  These  refined  models  included  a  (20-30-3)  architecture 
for  predicting  damping  in  the  first  three  modes,  a  (20-30-3)  model  for  predicting 
the  first  three  natural  frequencies,  a  (20-20-1)  model  for  predicting  the  static 
displacement,  and  a  (6-15-1)  model  for  predicting  the  failure  criterion.  Only  six 
input  neurons  were  used  in  the  input  layer  of  the  last  model,  as  only  6  design 
variables,  x^g  through  x^g,  representing  the  stacking  sequence  of  the  laminate  and 
its  ply  thickness,  were  required  to  determine  the  failure  criterion  of  the  laminate. 
Another  numerical  experiment  was  conducted  with  the  following  architectures, 
(20-30-20-3),  (20-30-30-3),  (20-20-20-1),  and  (6-20-20-1)  where  one  additional 
hidden  layer  was  added  in  the  previous  four  models.  These  four  models 
determined  3  dampings,  3  frequencies,  one  displacement  and  one  failure  criterion, 
respectively.  A  further  refinement  of  the  neural  network  models  which  compute 
damping  and  frequencies  were  also  built.  These  entailed  six  (20-30-20-1)  models, 
each  of  which  predicted  either  one  damping  or  one  frequency  in  a  specific 
vibrational  mode.  A  subscript  was  placed  after  the  parenthesis  representing  the 
specific  neural  network,  such  as  (20-30-30-8)ioo,  where  the  subscript,  100,  indicates 
the  total  number  of  input-output  training  sets  used  to  train  the  model. 
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The  finite  element  analysis  technique  [33]  was  used  to  compute  the 
structural  characteristics  of  the  three-dimensional  viscoelastically  damped  beam. 
The  beam  was  analyzed  as  having  plate  characteristics.  Plates  treated  with 
viscoelastic  layers  can  be  modeled  as  shown  in  Figure  6.6,  wherein  a  three  layer 
sandwich  is  modeled  by  two  layers  of  nodes.  Both  the  upper  and  lower  layers  can 
be  arbitrarily  laminated.  Two  hundred  and  twenty  randomly  generated  input 
vectors  were  created,  and  the  finite  element  analysis  approach  was  used  to 
generate  the  output  vectors  corresponding  to  these  inputs.  One  hundred  sets  of 
input-output  data  were  initially  used  to  train  selected  neural  models.  An 
additional  set  of  one  hundred  training  patterns  were  then  combined  with  the 
original  set  to  make  a  new  training  set  of  two  hundred.  The  two  hundred-set  data 
was  used  in  training  the  same  neural  models.  The  last  twenty  input-output 
mapping  data  sets  were  used  for  comparison  with  the  outputs  of  all  trained 
models  to  determine  the  accuracy  of  results  predicted  by  the  neural  networks.  All 
neural  network  training  was  done  on  a  Micro Vax  3300  computer.  The  CPU  time 
(minutes)  required  to  train  these  models  and  the  accuracies  of  the  mapping  by 
these  models  are  shown  in  Table  6.2. 

The  CPU  time  obviously  increased  when  training  was  performed  using  all 
200  input-output  patterns.  It  was  most  apparent  when  using  the  largest  model, 
(20-30-30-8),  where  fourteen  times  more  CPU  time  was  used  to  train  the  model 
with  200  training  pairs  than  with  100  pairs.  Smaller  models  were  observed  to 
require  less  CPU  time  in  completing  the  training,  and  provided  more  accurate 
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Figure  6.6  Damped  plate  finite  element  model. 
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mapping.  As  a  result  of  the  high  degree  of  nonlinearity  in  the  input-output 
relationships,  and  the  large  dimensionality  of  the  structural  system,  higher 
accuracies  were  difficult  to  obtain  without  a  larger  number  of  training  pairs. 
Better  results  of  training  of  these  models  will  be  obtained  by  choosing  narrower 
ranges  of  input  design  variable  values.  Although  the  average  error  of  these 
models,  when  calculating  the  values  of  frequencies  and  failure  status,  was  large, 
the  percentage  of  accurate  prediction  of  constraint  violations  was  quite 
acceptable.  Accurate  values  of  these  quantities  were  not  as  important  when 
designs  were  feasible. 

Best  models  for  predicting  damping  in  the  first  three  modes  were  the  3 
models  identified  as  (20-30-20- 1)200,  i"  which  each  determined  the  damping  of  a 
single  mode.  Similarly,  for  predicting  frequencies  the  models  which  had  the  least 
average  errors  were  from  the  three  (20-30-20- 1)200  models.  Best  predictions  of 
scaled  displacement  and  failure  status  were  models  (20-20-20- 1)200,  and  (6-20-20- 
1)200'  respectively. 

The  six  trained  neural  networks  mentioned  above  were  subsequently  used 
in  a  genetic  search  to  determine  approximations  of  damping,  frequencies, 
displacement  and  failure  status  for  the  structural  system.  The  genetic  search  was 
completed  after  50  generations  with  a  population  size  of  120.  The  crossover 
probability  P,  was  0.8,  and  the  mutation  probability  P„  was  0.01.  The  VE^^^ 
approach  described  in  Chapter  4  was  used  to  simultaneously  generate  optimal 
solutions  according  to  different  weighting  combinations  of  design  criteria.  The 
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final  designs  corresponding  to  the  four  combinations  of  weighting  parameters  are 
shown  in  Table  6.3,  and  the  configurations  of  the  damping  treatments  of  these 
final  designs  are  shown  in  Figure  6.7. 
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Table  6.3  Final  designs  of  genetic  search. 


damping. 

freq,  [Hz] 

X,,  through  Xn  fdeel 

w. 

damping. 

freq,  [Hz] 

X,o  [m] 

W3 

damping3 

freq3  [Hz] 

Xi9  [m] 
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weight  [N] 

disp  [m] 

X20  [m] 

failure 

0.5 

0.1529(0.1311)' 

4.0(4.3) 

20,  30,  30,  -75,  55 

0.0 

0.0834(0.0687) 

23.3(21.0) 

0.000160 

0.0 
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0.0006 

0.5 
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0.5 
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0.0 
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69.9(60.1) 
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0.5 
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0.2 

0.1506(0.1573) 
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0.2 
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0.000160 

0.2 

0.0414(0.0262) 
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0.4 

2.4154 
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-0.32(-0.15) 
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Figure  6.7  Damping  treatment  patterns  for  final  designs. 


CHAPTER  7 
CONCLUDING  REMARKS 


7.1  Conclusions 

The  present  study  has  examined  the  feasibility  of  using  genetic  search 
methods  in  optimization  problems  characterized  by  a  mix  of  continuous,  discrete, 
and  integer  design  variables.  Design  variable  representation  schemes  required  to 
facilitate  this  application  are  proposed  and  evaluated  for  their  performance  in 
representative  design  problems.  Genetic  search  methods  require  less 
computational  effort  than  the  traditional  branch-and-bound  schemes  in  problems 
where  the  number  of  integer  and  discrete  design  variables  is  large.  Additionally, 
the  methods  have  an  inherently  better  probability  of  locating  the  global  optimum 
than  the  gradient-based  mathematical  programming  methods. 

This  study  also  examines  the  issue  of  transforming  constrained  optimization 
problems  into  unconstrained  function  maximization,  as  required  in  the  genetic 
search  approach.  A  penalty  function  approach  was  used,  and  preliminary  findings 
underscore  the  need  for  a  very  careful  evaluation  of  the  objective  function 
magnitude  in  relation  to  that  of  the  penalty  term.  Further  numerical 
experimentation  on  large  scale  problems  is  necessary  to  build  upon  the  general 
observations  available  from  this  study.  It  was  found  that  caution  must  be 
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exercised  when  using  genetic  search  methods  to  solve  constrained  optimization 
problems. 

This  study  also  examines  the  use  of  genetic  algorithms  in  the  simultaneous 
generation  of  a  family  of  Pareto  optimal  designs.  Three  distinct  strategies  for 
implementing  such  a  procedure  were  proposed  and  examined  in  the  context  of 
structural  optimization  problems.  The  approach  was  shown  to  be  a  viable  option 
for  multicriterion  design  problems.  It  is  particularly  powerful  for  problems  with 
known  nonconvexities  in  the  design  space.  The  extensive  numerical  experience 
derived  in  this  implementation  was  considered  extremely  useful  for  the  further 
refinement  of  this  class  of  search  techniques. 

The  genetic  search  strategies  presented  in  this  study,  were  implemented  on 
a  multicriterion  optimization  problem  involving  the  optimal  design  of  a  two- 
dimensional,  viscoelastically  damped,  laminated  composite,  cantilever  beam.  In 
this  approach,  the  structural  damping  was  enhanced  by  a  combination  of 
redesigning  the  base  structure  and  by  application  of  a  viscoelastic  damping  tape 
on  the  surface  of  the  base  structure.  Results  of  the  optimization  problem  showed 
that  discrete  types  of  surface  treatment  of  viscoelastic  damping  tapes  were  favored 
over  a  continuous  type  of  treatment  when  the  damping  criteria  in  higher  modes 
were  heavily  weighted.  Discontinuous  damping  tape  on  the  beam  led  to  a 
noticeable  loss  in  transverse  stiffness  for  a  thin  beam.  Such  discrete  treatment  is 
not  preferred  when  a  tight  static  tip  displacement  constraint  is  in  effect. 
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A  neural-network-based  approximation  was  considered  in  this  work  to 
alleviate  the  computational  requirements  of  the  genetic  search  procedure.  The 
applicability  of  using  a  multilayer  perceptron  model  of  an  artificial  neural 
network,  to  map  the  input-output  relationships  between  the  design  variables  and 
the  output  response  of  a  three-dimensional,  viscoelastically  damped,  laminated 
composite  tapered  beam  was  examined.  Two  distinct  types  of  network  models 
were  considered.  The  first  category  (lumped  network)  was  one  where  the  input 
patterns  were  mapped  into  all  of  the  required  output  quantities.  The  second 
approach  (distributed  network)  was  based  on  creating  several  subnetworks,  where 
each  subnetwork  mapped  the  input  quantities  into  a  smaller  subset  of  the  output 
quantities.  The  training  of  distributed  models  required  less  computational  effort 
and  resulted  in  higher  mapping  accuracies  than  the  lumped  network  models.  This 
was  even  more  significant  as  the  number  of  training  patterns  needed  for  the 
complex  structural  system  were  increased.  Furthermore,  a  group  of  distributed 
models  with  a  relatively  small  number  of  output  elements  is  advantageous  when 
the  structural  system  is  re-optimized  with  a  different  set  of  objective  functions  and 
constraints.  In  the  evem  that  an  optimal  design  is  to  be  generated  with  additional 
constraints,  a  small  model  must  be  trained  which  contains  this  constraim  as  an 
output  node.  Additionally,  if  the  definition  of  one  constraim  is  changed,  it  only 
requires  that  the  small  network  be  retrained.  In  the  case  of  generating  a  design 
with  fewer  constraints,  the  optimization  can  proceed  by  simply  discarding  the 
models  which  map  the  input  vector  into  the  deleted  output  elements.  A  lengthy 
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training  process  of  a  new  neural  network  modeling  an  entire  input-output 
relationship  of  the  structural  system  would  otherwise  be  necessary.  It  is  also 
worthwhile  to  indicate  that  the  training  times  can  be  reduced  further  by  training 
distributed  models  in  parallel  if  such  a  computational  architecture  is  available. 

7.2  Recommendations  for  Future  Work 
The  exterior  penalty  function  approach  was  exclusively  used  in  the  genetic 
search  algorithms  to  transform  constrained  optimization  problems  into 
unconstrained  optimization  problems  in  this  study.  Different  approaches  should  be 
examined  with  the  objective  of  seeking  a  possibly  better  approach  to  deal  with 
constrained  optimization  problems.  These  include  different  types  of  control 
schemes  to  limit  the  magnitude  of  the  penalty  term.  In  the  present  work,  it  was 
found  that  without  such  explicit  limits,  the  convergence  rate  was  significantly 
'\  lowered.  ■  " "  ' 

In  the  context  of  multicriterion  optimization,  the  use  of  a  population 
average  value  of  a  specific  criterion  for  a  given  generation  may  be  replaced  by 
other  schemes.  Caution  is  necessary  in  the  sharing  approach  for  multicriterion 
optimization,  as  it  is  analogous  to  the  problem  of  biological  in-breeding.  Sources 
of  diversity  in  the  population  must  be  introduced.  Mating  restrictions,  other  than 
those  based  entirely  on  the  weighting  coefficient  combinations,  will  be  important 
to  include  in  future  work.  Since  the  population  of  designs  in  each  specific 
weighting  group  is  based  on  the  penalized  fitness  of  the  entire  population,  a  better 
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average  fitness  for  a  specific  weighting  group  will  eventually  result  in  that  group 
attracting  more  members  in  successive  generations.  The  fitness  or  the  overall 
utility  function  has  to  be  carefully  scaled  so  as  to  distribute  designs  evenly  over 
each  weighting  combination,  and  for  each  weighting  group  to  possess  roughly 
similar  exploitative  abilities  in  the  genetic  search. 

_    Due  to  the  possible  small  size  of  each  subpopulation  in  the  VEj^c  and 
^Ee,.c  approaches,  genetic  drift  might  result,  with  associated  undesirable 
characteristics  such  as  convergence  to  a  sub-optimal  design.  A  possible  approach 
to  prevent  such  a  premature  convergence  is  to  introduce  a  more  moderate  sized 
population  for  genetic  evolution.  The  underlying  assumption  here  is  that  a  larger 
population  is  likely  to  be  more  diverse,  and  hence  offers  a  higher  possibility  of 
locating  the  true  optimum.  This  can  be  achieved  in  a  strategy  where  the 
crossover  transformation  is  a  two-stage  process.  In  the  first  stage,  a  crossover  is 
performed  between  a  subset  of  members  within  the  same  subpopulation  (akin  to 
inbreeding).  The  remaining  members  from  each  subset  form  a  larger  gene  pool 
where  a  second  stage  of  crossover  is  conducted.  In  this  second  stage,  a  mating 
restriction  must  be  applied  on  the  variable  denoting  the  weighting  combination 
associated  with  a  design.  This  would  prevent  crossovers  between  designs 
corresponding  to  radically  different  weight  combinations. 

When  a  neural  network  is  used  to  approximate  the  input-output 
relationship  of  a  complex  structural  system,  a  relatively  large  number  of  training 
pairs  may  be  necessary.  In  this  situation,  distributed  neural  network  models  with  a 
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smaller  number  of  elements  in  the  output  layer  were  shown  to  be  effective  in 
saving  computational  effort.  Additional  experiments,  where  a  group  of  distributed 
models  is  chosen  to  replace  a  large  model  are  necessary  to  establish  more 
fundamental  understanding  of  the  method.  Another  approach  which  employs 
clustering  techniques  [98]  to  narrow  the  range  of  input  vectors  may  be  used  for 
complex  structural  systems  in  future  work.  A  clustering  approach  could  also  be 
used,  in  future  work,  in  conjunction  with  the  refinement  approach  presented  in 
this  work.  Further  research  is  required  to  investigate  approximation  schemes 
other  than  the  neural  network  based  approximation  used  in  this  work.  These 
methods  may  include  linear  approximation  based  on  first  order  sensitivity,  and  the 
accumulated  approximation  techniques  proposed  by  Rasmussen  [103]. 
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